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This  dissertation  Is  the  result  of  my  application  of  quaslllnear 
plasma  theory  to  the  high  altitude  nuclear  explosion  debris-air  coupling 
problem.  Although  I have  written  the  dissertation  with  this  problem  In 
mind,  the  plasma  theory  and  numerical  techniques  are  applicable  to  a wide 
range  of  physical  problems  In  geophysics  and  In  controlled  thermonuclear 
reactor  research.  I have  written  the  dissertation  with  both  the  reader 
familiar  with  quaslllnear  theory  and  the  reader  with  only  a limited 
knowledge  of  plasma  physics  In  mind.  The  experienced  reader  can  omit 
the  physical  discussion  sections  in  Chapter  II  and  all  of  the  derivations 
In  the  Appendices  without  loss  In  continuity.  These  discussions  are 
primarily  for  the  Inexperienced  reader  wishing  to  learn  quaslllnear 
theory.  I have  assumed  that  this  reader,  who  wishes  to  understand  all 
of  the  derivations,  has  some  basic  knowledge  of  plasma  physics  including 
Vlasov  theory.  If  not,  the  references  given  throughout  the  dissertation 
contain  the  necessary  Information.  For  the  reader  Interested  In  the 
computer  program,  I have  Included  an  outline  of  the  entire  program  In 
Appendix  G. 
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Notation 


atomic  weight  (atomic  mass  units) 
magnetic  field  (gauss) 

O O 

diffusion  coefficient  (cm  /sec  ) 
electric  field  (statvolt/cm) 

spectral  wave  energy  (ergs  in  3D  or  ergs/cm  In  2D) 

p 

current  density  (statcoulomb/cm  -sec) 

p 

pressure  (dynes/cm  ) 
temperature  (degrees  Kelvin) 
volume  (cm  ) 

computer  velocity  mesh  translation  constant  (cm/sec) 

Initial  relative  velocity  between  the  debris  and  air  ion 
beams  (cm/sec) 

species  average  velocity  (cm/sec) 

3 

energy  (ergs)  or  energy  density  (ergs/cm  ) 
plasma  dispersion  function  - 

7m  M oc 

speed  of  light  (cm/sec)  * 3.00(10*®) 
thermal  speed  of  species  <*  (cm/sec) 

thermal  speed  of  species  * In  wave  propagation  direction 
(cm/sec) 

elementary  charge  (statcoulomb)  * 4.8(10"*°) 
distribution  function  (sec°/cm°  In  3D  or  sec2/cm^  in  2D) 

yr 

wave  propagation  vector  (cm-*)  * 2TK/wavelength 
mass  (gm) 

number  density  (cm  v) 
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V Uz 

vxo*  vzo 


S Le 

J*\.\ 


number  of  elementary  charges  (q  = 1 means  singly  charged) 
time  (sec) 

computer  velocity  mesh  coordinates 

computer  velocity  mesh  scale  constants  (cm/sec) 

particle  velocity  (cm/sec) 


ion  acoustic  speed  (cm/sec) 
Alfve'n  wave  speed  (cm/sec)  ! 


_ /Ye  XT, 


vph  wave  phase  velocity  (cm/sec)  =-^-^4 
~x  particle  location  (cm) 

x,  y,  z particle  location  (cm) 

<*  denotes  particle  species  (debris,  air,  or  electrons) 


electron  beta  = 


//f 


ciemuM  ucl a - 

ratio  of  specific  heat  at  constant  pressure  to  that  at 
constant  volume 

plasma  dielectric  of  species  °< 

angle  between  the  streaming  direction  and  the  wave  propaga- 
tion direction  (degrees) 


Boltzmann  constant  (erg/degree  Kelvin)  = 1.38(10”  ) 

I J ^ x 

electron  gyro  radius  vs  wavelength  parameter  = 

Dc-bye  length  for  species  <<  (cm)  = 


anomalous  collision  frequency  (sec  ) = — 


pi  = 3.14 


^a(  VA  ~ vt> 


average  electron  gyro  (Larmor)  radius  (cm)  = y 

\/  • 

Ion  gyro  radius  (cm)  = — . 

time  normalization  constant  (sec)  = 

angle  between  velocity  point  of  Interest  and  wave 
propagation  direction  (degrees) 

complex  wave  frequency  (radians/sec) 
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_n 

« a.  / A l 1 B 1 

gyro  frequency  (sec  A)  ■ 4 /vn  e. 

■ 

plasma  frequency  (sec  ^ 

lower  hybrid  frequency  (cm  ) ■ + _Ji\ 

~LH 

0 subscript 

denotes  zero  order  quantity 

1 subscript 

denotes  first  order  quantity 

1 subscript 

imaginary  part  or  Ion  species 

j- 

r subscript 

D,  A,  or  e 

real  part 

• 

■ 

subscript 

species:  debris,  air,  or  electrons 

1 i 

i i 

c(  subscript 

species 
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Abstract 

Improved  calculations  of  the  anomalous  transport  coefficients  due 
to  the  Ion-Ion  two-stream  plasma  Instability  are  presented.  The  coeffi- 
cients describe  the  transfer  of  momentum  and  energy  between  the  Ion 
species  by  the  electric  fields  generated  by  unstable  waves.  The  calcula- 
tions use  electrostatic  quaslllnear  theory  In  two  dimensions  to  trace 
the  time  evolutions  of  the  two  ion  species.  The  evolutions  are  continued 
until  the  transport  coefficients  reach  their  maximum  values.  The  evolu- 
tions are  done  for  homogeneous  plasmas  with  warm  Ion  beams  and  a single 
hot  electron  background  distribution.  The  electrons  may  or  may  not  be 
magnetized  by  a uniform  magnetic  field  perpendicular  to  the  Ion  streaming 
motion.  Because  the  calculations  are  two  dimensional,  the  unstable  ion- 
ion  waves  do  not  necessarily  propagate  In  the  streaming  direction.  The 
presentation  includes  derivations  of  quaslllnear  theory  from  Vlasov's 
and  Poisson's  equations  and  several  physical  discussions  regarding  the 
quaslllnear  equations.  The  analytical  forms  of  the  anomalous  transport 
coefficients  are  also  derived.  A detailed  discussion  of  the  computer 
program,  IONION,  which  was  written  to  numerically  solve  the  quaslllnear 
equations  Is  included.  Several  test  evolutions  are  done  for  plasmas 
which  could  characterize  the  early  time  nuclear  debris  cloud  surrounding 
a nuclear  detonation  In  the  terrestrial  exosphere.  The  two  Ion  species 
represent  the  nuclear  debris  Ions  In  the  expanding  cloud  and  the  Ionized 
ambient  air  being  overrun  by  the  debris.  The  results  of  the  evolutions 
Indicate  that  a fair  amount  of  momentum  transfer  from  the  debris  to  the 
air  occurs  when  the  relative  velocity  between  them  Is  less  than  the 
characteristic  wave  phase  velocity  of  the  plasma,  the  Ion  acoustic  speed 
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or  the  Alfven  speed  (respectively  for  unmagnetized  or  magnetized  elec 
trons).  On  the  other  hand,  little  transfer  occurs  when  the  relative 
velocity  Is  higher  than  the  characteristic  speed. 
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ANOMALOUS  TRANSPORT  COEFFICIENTS  DUE  TO  THE 
ION-ION  TV/O- STREAM  PLASMA  MICROINSTABILITY 
I.  Introduction 

Exospheric  nuclear  bursts  Inject  large  numbers  of  radioactive 
fission  debris  nuclei  into  the  terrestrial  magnetosphere.  These  radio- 
active nuclei  decay  emitting  copious  numbers  of  energetic  electrons 
which  often  are  trapped  In  the  magnetosphere.  These  energetic  electrons 
(beta  particles)  supplement  the  naturally  occurring  Van  Allen  radiation 
belts  greatly  enhancing  the  flux  (electrons/cm  /sec)  of  high  energy 
electrons.  (A  general  discussion  of  the  magnetosphere  and  trapped  radia- 
tion Is  In  Ref  1.)  Because  the  beta  electrons  typically  have  energies 
of  the  order  of  1 MeV,  satellites  passing  through  the  enhanced  belts  can 
be  damaged.  This  damage  problem  is  particularly  severe  to  the  large 
scale  Integrated  circuitry  being  introduced  on  current  and  future 
satellites  as  these  components  often  have  a relatively  low  radiation 
damage  threshold.  The  problem  Is  of  great  Interest  to  the  agencies  who 
manage  satellite  systems,  and  much  effort  has  been  expended  to  Improve 
the  beta  radiation  prediction  capability. 

The  topic  of  this  dissertation  Is  an  outgrowth  of  the  following 
facet  of  the  beta  electron  Injection  process:  the  possible  mixing  of 
the  nuclear  debris  with  exospheric  air  Ions  by  the  ion-ion  Instability. 
This  Instability  could  occur  during  the  initial  outward  expansion  of  the 
debris.  This  expansion  lasts  about  a second  after  the  detonation  of  the 
nuclear  device.  If  the  Instability  occurs  and  successfully  causes  the 
mixing,  the  air  would  accelerate  to  the  debris  expansion  velocity. 

This  acceleration  of  the  air  Is  called  debris-air  coupling.  The  coupling 
problem  Is  discussed  in  more  detail  In  the  next  section  of  this  chapter. 
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In  the  sections  of  this  chapter  following  that  discussion,  the  disserta- 
tion calculation  of  the  ion-ion  Instability  Is  outlined.  Although  the 
dissertation  addresses  only  the  Ion-Ion  instability  in  the  context  of 
debris-air  coupling,  the  calculatlve  procedure  Is  applicable  to  a wide 
range  of  physical  problems  in  the  fields  of  geophysics  and  controlled 
thermonuclear  reactor  research. 

Debris-Air  Coupling  Problem 

The  SPECTER/SAFER  computer  program,  which  has  been  developed  jointly 
by  the  Air  Force  Weapons  Laboratory  (AFWL)  and  the  Lockheed  Palo  Alto 
Research  Laboratory,  calculates  the  accumulated  radiation  fluence  In 
electrons  per  square  centimeter  received  by  a satellite  from  a given 
exospheric  nuclear  burst.  (The  base  of  the  exosphere  is  considered  to 
be  at  an  altitude  of  100  km.)  The  program  first  calculates  the  initial 
expansion  of  the  nuclear  debris  as  It  pushes  against  the  residual  atmo- 
sphere and  the  magnetic  field  of  the  earth  to  form  a debris-air  bubble. 
The  cross  section  of  a typical  bubble  which  surrounds  the  burst  point  Is 
shown  in  Fig.  1.  The  program  then  calculates  the  extension  of  the  Ion- 
ized debris  along  the  magnetic  field  lines  affected  by  the  bubble.  This 
extension  proceeds  until  the  majority  of  the  debris  following  the  field 
lines  ends  up  in  the  upper  atmosphere  of  the  earth  at  the  conjugate 
regions  shown  in  the  figure.  The  program  calculates  the  beta  decay  of 
the  debris  as  it  moves  and  then  calculates  the  formation  of  the  radiation 
belts.  Finally,  the  program  computes  the  fluences  received  by  satellites 
moving  through  the  belts. 

One  of  the  factors  causing  uncertainties  in  the  final  fluences 
calculated  by  the  program  Is  the  uncertain  fraction  of  the  ambient  air 
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Fig.  1.  The  Formation  of  a Debris-Air  Bubble  In  the  Exosphere, 
during  the  Initial  bubble  formation  Is  the  result  of  a transfer  of 
momentum  and  energy  from  the  expanding  debris  ions  to  the  air  atoms  which 
have  been  Ionized  by  x-rays  and  ultraviolet  radiation  from  the  bomb. 

If  the  amount  transferred  Is  small,  the  debris  expands  outward 
Impeded  only  by  the  magnetic  field  forming  a relatively  large  bubble.  If 
the  reverse  holds,  most  of  the  air  Is  accelerated  to  the  expansion 
velocity.  This  accumulation  of  moving  mass  rapidly  slows  the  outward 
expansion  so  the  bubble  Is  smaller.  The  resulting  radiation  belts  are 
then  smaller  and  more  Intense  significantly  affecting  the  fluences 
received  by  orbiting  satellites.  For  example,  a sample  calculation  has 
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shown  that  the  radius  of  the  bubble  from  a one  klloton  detonation  at  400 
kilometers  Is  cut  In  half  when  the  fraction  of  air  coupled  to  the  debris 
Is  Increased  from  zero  to  unity  (Ref  2). 

The  actual  fraction  of  debris-air  coupling  Is  calculated  by  KLYSMA, 
a large  hydrodynamic  computer  program  written  by  the  Naval  Research 
Laboratory  (NRL).  (Examples  on  the  nature  of  hydrodynamic  theory  and 
some  references  on  the  subject  are  In  Chapter  II  of  the  dissertation.) 

The  coupling  cannot  be  due  to  conventional  binary  collisions  In  which 
the  debris  and  air  Ions  simply  collide  with  each  other.  At  the  air  Ion 
densities  characteristic  of  even  a fully  Ionized  exosphere,  the  average 
distance  an  Individual  debris  ion  travels  between  binary  collisions  Is 
of  the  order  of  hundreds  of  kilometers.  This  distance  is  larger  than  a 
typical  debris-air  bubble.  However,  the  coupling  can  be  due  to  collective 
Interactions  between  electromagnetic  fields  and  the  plasma  species: 
debris,  air,  and  electrons.  In  collective  interactions,  large  numbers  of 
the  plasma  particles  collectively  participate  to  transfer  momentum  and 
energy  among  them  via  the  fields.  The  program  KLYSMA  Includes  several 
of  these  Interactions  by  incorporating  anomalous  transport  coefficients 
Into  Its  hydrodynamic  equations.  The  term  "anomalous"  means  the  trans- 
port coefficients,  which  appear  In  the  hydrodynamic  equations  in  the  same 
manner  as  those  due  to  ordinary  binary  collisions,  actually  refer  to 
collective  Interactions. 

One  of  the  most  Important  of  the  collective  Interactions  causing 
coupling  Is  the  Ion-Ion  two-stream  plasma  mlcroinstabllity  (Ref  3:  1-2). 

In  a plasma  Instability,  available  energy  from  a nonequlllbrlum  feature 
Is  explosively  converted  Into  rapidly  growing  electromagnetic  waves. 

These  waves  then  alter  the  plasma  tending  to  drive  It  toward  an  equlllb- 
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rlum  configuration.  The  energy  source  for  the  Ion-Ion  two-stream  Insta- 
bility Is  the  relative  velocity  between  the  debris  and  air  Ions.  The 
two  Ion  species  form  the  streams  or  beams  of  this  Instability,  and  the 
growing  waves  tend  to  reduce  their  relative  velocity.  The  term  "micro- 
instability"  means  the  external  appearance  of  the  plasma  Is  unchanged 
by  this  Instability.  Rather,  the  Instability  alters  the  velocity  dis- 
tributions of  the  Ions.  (A  velocity  distribution  describes  the  probabll- 

Ilty  that  a plasma  particle  Is  traveling  at  a given  velocity.  Typically, 
the  probability  distribution  Is  In  a gausslan  shape— the  Maxwell -Boltzman 
distribution.  The  average  velocity  of  the  gausslan  Is  the  beam  velocity, 
and  the  spread  around  the  mean  Is  the  thermal  velocity  which  Is  related 
to  the  temperature.  Each  plasma  species  has  Its  own  distribution,  and 
the  Integral  of  Its  distribution  over  all  possible  velocities  equals  the 
species  density  In  particles  per  unit  volume.) 

Dissertation  Calculation 

The  anomalous  transport  coefficients  from  the  Ion-Ion  Instability 
have  been  Included  In  KLYSMA.  The  coefficients  are  obtained  from  one- 
dimensional  calculations  which  assume  all  of  the  growing  electromagnetic 
fields  generated  by  the  Instability  propagate  In  the  radial  direction 
(the  local  expansion  direction).  For  reasons  which  are  discussed  later, 
these  one-dimensional  calculations  are  inadequate  for  large  expansion 
♦ velocities. 

Thus,  the  objective  of  this  dissertation  is  to  calculate  improved 
anomalous  transport  coefficients  due  to  the  Ion-Ion  two-stream  plasma 
mlcrolnstablllty.  The  plasma  parameters  used  In  the  calculations  such 
as  the  species  densities,  magnetic  field  strengths,  temperatures,  and 

w 

relative  velocities  are  chosen  to  be  representative  of  those  which  could 
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be  found  In  the  solution  of  a debris-air  coupling  problem.  The  calcula- 
tions are  done  simulating  both  large  and  small  expansion  velocities.  In 
the  next  few  sections,  the  general  nature  of  the  dissertation  calculations 
Is  outlined.  First,  the  method  used  to  obtain  the  anomalous  coefficients 
Is  described.  Next,  the  type  of  plasma  and  the  assumptions  used  are 
outlined.  And  finally,  how  the  dissertation  advances  the  state  of  the 
art  Is  explained. 

Calculatlve  Method.  The  calculatlve  method  Is  to  solve  for  wave- 
particle  Interactions  the  set  of  electrostatic  quasll Inear  equations 
described  In  the  next  chapter.  In  quasi! inear  theory,  a spectrum  of  weak, 
unstable  waves  Is  assumed  to  slowly  modify  the  spatially  averaged  velocity 
distributions  of  the  plasma  species.  As  the  averaged  distributions 
change,  the  oscillation  frequency  and  wave  amplitude  growth  rate  of  each 
of  the  unstable  wave  modes  in  the  spectrum  slowly  change.  Eventually, 
the  modifications  saturate  when  the  waves  can  no  longer  change  the 
distributions. 

The  analytical  set  of  equations  is  solved  numerically  on  a Control 
Data  Corporation  (CDC)  7600  computer  using  a computer  program  expressly 
written  for  this  task.  A particle  simulation  calculation  is  not  done 
because  the  analytical  approach  saves  computer  time  and  gives  better 
physical  Insight.  The  computer  solution  follows  the  changes  of  the 
species  distributions  with  time.  This  following  of  the  changes  In  time 
Is  called  the  evolution  of  the  plasma.  Periodically  during  the  evolution, 
the  wave  mode  frequencies  and  growth  rates  are  updated  by  solutions  of 
the  dispersion  relation,  one  of  the  equations  described  in  the  next 
chapter.  The  calculation  Is  called  quaslllnear  because  linear  plasma 
theory,  the  solutions  of  the  dispersion  relation,  is  used  not  only  for 
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the  Initial  species  distributions  but  also  for  distributions  which  have 
been  evolved  In  time. 

This  dissertation  Includes  the  evolution  of  the  species  distributions 
for  plasmas  which  generate  unstable  Ion-Ion  waves.  The  distributions  are 
evolved  In  time  until  the  growing  transport  coefficients,  which  are 
related  to  the  wave  energies,  peak  and  then  begin  to  decline  Indicating 
the  approach  to  quaslllnear  saturation.  Saturation  of  the  ion-ion  insta- 
bility can  occur  In  the  debris-air  coupling  problem  because  the  typical 
Instability  growth  time  of  about  one  microsecond  is  much  shorter  than 
the  typical  bubble  expansion  time.  The  transport  coefficients  return  to 
zero  at  saturation,  but  quasilinear  theory  is  then  no  longer  valid.  The 
dissertation  calculations,  therefore,  give  upper  limits  to  the  coeffi- 
cients to  help  determine  whether  debris-air  coupling  is  possible  and  to 
help  estimate  the  coupling  rate.  The  evolutions  are  necessary  to 
determine  the  maximum  values  of  the  transport  coefficients  because  these 
values  cannot  be  accurately  estimated  from  the  Initial  growth  rates  of 
the  unstable  waves.  Two  different  plasmas  generating  waves  with  the 
same  growth  rates  can  couple  quite  differently.  An  example  of  this  is 
in  Chapter  IV. 

Plasma  Type.  The  transport  coefficients  are  found  for  the  following 
type  of  plasma:  (1)  The  species  distributions  are  initially  spatially 
homogeneous  so  the  calculation  represents  what  is  taking  place  in  a 
small  volume  in  the  mixing  region  near  the  outer  edge  of  the  debris 
bubble.  (See  Fig.  1.)  The  few-kllometer- thick  mixing  region  is  much 
smaller  than  the  bubble  size  and  Is  where  the  debris-air  mixing  takes 
place.  The  homogeneous  approximation  Is  relatively  good  when  the  un- 
stable mode  wavelengths  and  the  coupling  length  (the  distance  the  debris 


AFIT/DS/PH/78-1 


"7 

\ I 

I 

moves  during  the  evolution  to  saturation)  are  much  less  than  the  mixing 
region  size.  (2)  There  are  no  ambient  electric  fields.  (3)  A uniform 
magnetic  field  is  In  the  z direction  which  is  perpendicular  to  the  debris 

streaming  in  the  x direction.  If  the  debris  expansion  velocity  is 

f * 

approximately  equal  to  or  greater  than  the  Alfven  speed,  the  magnetic 
field  of  the  earth  Is  nearly  perpendicular  to  the  local  expansion 
direction  nearly  everywhere  in  the  mixing  region  which  forms  the  edge  of 
the  bubble.  This  occurs  because  the  Alfven  speed  is  the  speed  a low 
frequency  disturbance  can  propagate  In  a magnetized  plasma.  (This  speed 
Is  discussed  further  in  Chapter  II.)  The  bubble  expansion  Is  such  a low 

' • 

frequency  disturbance,  and  as  a result,  the  terrestrial  field  often  piles 
up  on  the  outside  of  the  bubble  as  shown  in  Fig.  1.  (The  field  cannot 
easily  penetrate  the  interior  of  the  bubble  because  It  Is  filled  with  a 
hot,  nearly  perfectly  conducting  plasma.)  Therefore,  the  chosen  plasma 
magnetic  field  configuration  can  reasonably  represent  the  compressed 
field  of  the  earth  In  a small  volume  in  the  mixing  region.  (4)  The 
debris  and  air  Ions  are  assumed  to  be  initially  in  separate  Maxwellian 
velocity  distributions  with  the  relative  velocity  between  them  being  the 
bubble  expansion  velocity.  These  distributions  are  "warm,"  meaning  that 
the  Ion  thermal  velocities  may  be  comparable  to  the  relative  velocity. 
Also,  the  Ion  distributions  may  each  have  a separate  temperature  in  the 
x and  z directions.  The  densities,  atomic  weights,  and  ionization  levels 
of  the  Ion  species  are  not  restricted  although  the  densities  will  usually 
be  comparable  because  experience  has  shown  that  most  of  the  coupling 
occurs  at  that  time  (Ref  4).  (5)  All  of  the  electrons  are  assumed  to  be 

In  a single  background  Maxwellian  distribution  possibly  with  different 
x and  z direction  temperatures.  These  electrons  from  both  ion  species 
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are  assumed  to  have  been  mixed  together  by  some  faster  acting  process 
such  as  an  electron-ion  instability.  Their  distribution  thermal  veloci- 
ties will  be  larger  than  the  ion  relative  velocities  in  agreement  with  the 
results  of  electron-ion  instability  computer  simulations  (Refs  5 and  6). 
(6)  The  dissertation  calculation  assumes  that  each  ion  species  is  fully 
ionized  at  a given  ionization  level.  No  neutral  particles  are  included. 
These  approximations  are  reasonable  because  most  of  the  ions  of  each 
species  in  a hot  plasma  tend  to  be  at  the  same  ionization  level  and 
because  collisions  are  Ignored  so  neutral  particles  cannot  influence  the 
ions. 

Calculation  Assumptions.  Several  basic  assumptions  are  made  in  the 
solution  of  the  quaslllnear  equations.  First,  the  equations  are  assumed 
to  be  valid,  meaning  that  the  species  distributions  change  relatively 
slowly  because  the  wave  electromagnetic  fields  are  relatively  weak.  This 
condition  permits  the  use  of  perturbation  theory  because  the  changes  in 
the  distributions  during  their  evolutions  are  always  perturbations  to  a 
known  zero  order,  time  Independent  solution.  (See  Ref  7:  66-67  for  a 
typical  example  of  perturbation  theory.)  Quasilinear  theory  neglects 
small,  high  order  effects  such  as  the  coupling  of  energy  from  one  wave 
mode  into  another.  Each  wave  mode  is  assumed  to  remain  independent  from 
all  of  the  rest. 

The  calculation  is  done  In  two  dimensions.  The  velocity  space 
distributions  of  the  ions  are  evolved  in  the  vx-vz  plane  interacting 
with  plane  waves  having  propagation  directions  in  the  kx-kz  plane.  Slab 
geometry,  not  cylindrical  symmetry,  has  been  assumed  for  the  ion  distri- 
butions and  the  waves  so  their  coordinates  are  strictly  Cartesian.  This 
means  the  shape  of  the  ion  distributions  and  the  nature  of  the  wave 
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spectrum  are  unchanged  by  a translation  In  the  y direction.  The  calcula- 
tion would  have  been  one  dimensional  if  the  waves  propagated  only  in  the 
direction  and  interacted  with  distributions  which  were  functions  of 
vx  only.  The  reason  why  the  calculation  has  to  be  done  in  two  dimensions 
is  that  for  large  relative  velocities  the  unstable  wave  modes  propagate 
nearly  perpendicular  to  the  ion  streaming  direction. 

The  electrons  are  "adiabatic,"  meaning  their  distribution  does  not 
change  with  time.  The  Initial  type  of  plasma  is  chosen  so  the  ion-ion 
Instability  dominates  over  electron-ion  instabilities  such  as  the  ion 
acoustic  or  the  modified  two-stream  Instability.  The  ion-ion  instability 
does  not  affect  the  electron  distribution  so  it  always  is  the  original 
Maxwellian. 

The  two  dimensionality  of  the  calculation  is  relaxed  for  the  elec- 
trons when  the  magnetic  field  significantly  affects  their  motions.  In 

these  magnetized  cases,  their  motions  around  the  field  in  the  v -v 

x y 

plane  are  averaged  assuming  cylindrical  symmetry.  This  process  is  valid 
when  the  wave  frequencies  are  much  less  than  the  electron  gyro  frequency. 
(A  gyro  frequency,  which  Is  defined  in  the  next  chapter,  is  the  rate  a 
charged  particle  orbits  around  a magnetic  field  line.  See  also  Ref  7: 
153-156.)  Because  the  electrons  are  adiabatic,  the  only  change  to  the 
quasi linear  equations  when  the  field  becomes  significant  is  to  the  elec- 
tron dielectric  term  In  the  wave  dispersion  relation. 

The  calculation  is  assumed  to  be  electrostatic,  meaning  the  magnetic 
field  of  each  wave  vanishes  because  the  curl  of  Its  electric  field  is 
zero.  The  electric  field  of  an  electrostatic  wave  Is  always  in  the  wave 
propagation  direction.  For  the  magnetized  cases,  a test  criterion  is 
evaluated  to  see  if  the  neglect  of  the  electromagnetic  terms  in  the 
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dispersion  relation  introduces  a significant  amount  of  error.  Electro- 
static wave  modes  often  grow  much  faster  than  electromagnetic  modes  so 
they  can  be  expected  to  dominate  the  evolution  (Ref  3:  29). 

The  ions  are  always  assumed  to  be  "unmagnetized,"  meaning  they  do 
not  see  the  magnetic  field.  This  approximation  is  a good  one  as  long  as 
the  wave  frequency  is  much  greater  than  the  ion  gyro  frequencies,  and  the 
ion  gyro  radii  are  much  greater  than  the  mode  wavelengths  and  the 
coupling  length. 

Relativistic  effects  are  ignored.  The  relative  velocities  of 
interest  are  at  least  two  orders  of  magnitude  less  than  the  velocity  of 
light. 

Advancement  of  the  State  of  the  Art.  The  complete  two-dimensional 
quasillnear  evolutions  which  are  done  in  this  dissertation  advance  the 
state  of  the  art.  Calculations  for  the  ion-ion  instability  have  been 
done  for  two  cold,  equally  dense  beams  of  the  same  ion  species  (Refs  8 
and  9).  The  calculations  in  Reference  8 are  done  in  one  dimension  using 
hydrodynamic  theory  which  ignores  the  effects  of  the  species  temperatures. 
Stability  is  predicted  for  large  relative  velocities.  The  computer 
simulation  calculations  in  Reference  9 are  done  for  a magnetic  field  free 
plasma  and  for  a ion  to  electron  mass  ratio  of  25.  The  Landau  damping 
and/or  growth  effects  (which  are  explained  in  the  next  chapter)  on  the 
waves  due  to  the  nonzero  thermal  velocities  of  the  ion  distributions  are 
not  adequately  included  in  either  of  these  results.  Conversely,  the 
dissertation  calculations  emphasize  these  effects.  Also,  the  dissertation 
calculations  can  be  done  for  cases  where  the  two  ion  beams  have  different 
densities,  types  of  particles,  temperatures,  and  ionization  levels. 

Cases  where  the  electrons  are  magnetized  or  unmagnetized  can  be  evolved 
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In  two  dimensions.  This  two  dimensionality  of  the  dissertation  calcula- 
tions permits  the  use  of  unstable  wave  modes  which  propagate  In  directions 
other  than  the  streaming  direction  and  permits  the  following  of  the  time 
evolutions  of  the  species  distributions  In  two  dimensions. 

Interestingly,  before  the  quaslllnear  evolutions  could  be  done, 
linear  theory  (the  use  of  the  dispersion  relation  for  the  Initial  distri- 
butions) had  to  be  developed  In  sufficient  generality  to  be  of  use  In  the 
quaslllnear  evolutions.  A multl-dimenslonal  magnetic- field- free  disper- 
sion relation  solution  for  warm  Ion  beams  has  been  done  (Ref  10). 

Certainly,  no  large  scale  numerical  solution  of  the  quaslllnear 
equations  has  been  done  to  Include  as  many  effects  as  are  Included  In  the 
dissertation  calculations.  The  program  written  to  do  these  calculations 
was  a major  undertaking,  and  none  like  It  exists.  These  calculations 
evolve  particles  both  resonant  and  nonresonant  with  the  waves  (Resonant 
particles  have  velocities  approximately  equal  to  the  wave  phase  velocities 
as  will  be  explained  later.).  Include  the  temperatures  of  the  Ion  beams. 
Include  the  electron  temperature,  account  for  two  dimensionality,  and 
Include  an  ambient  magnetic  field.  The  availability  of  a computer  with 
the  speed  of  the  CDC  7600  so  that  the  millions  of  necessary  Integrations 
could  be  done  in  a reasonable  amount  of  time  was  extremely  helpful  In 
allowing  this  degree  of  completeness. 

The  above  discussion  regarding  the  advancement  of  the  state  of  the 
art  completes  the  Introduction  of  the  topic  of  this  dissertation.  First, 
the  background  motivation  of  this  topic,  the  exospheric  nuclear  detona- 
tion debris-air  coupling  problem,  was  presented.  Next,  the  dissertation 
objective  and  the  calculative  methods  used  to  reach  that  objective  were 
described.  Then,  the  type  of  plasma  and  the  calculation  assumptions 
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were  presented.  And  finally,  why  this  dissertation  Is  a significant 
advance  In  the  field  of  plasma  physics  was  explained. 

The  remainder  of  the  dissertation  Is  as  follows:  In  the  next 
chapter,  the  equations  of  quaslllnear  theory  are  presented,  and  their 
physical  significances  are  explained.  In  Chapter  III,  the  solution 
procedure  Is  described.  In  Chapter  IV,  the  results  of  the  dissertation 
calculations  are  presented,  and  their  physical  significances  are 
explained.  Finally,  In  the  last  chapter,  the  results  are  summarized  and 
recommendations  for  further  study  are  made. 
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II.  The  Equations  of  Quaslllnear  Theory 


In  this  chapter,  the  equations  for  electrostatic  quaslllnear  wave- 
particle  Interactions  are  presented.  Most  of  the  algebraic  details  fn 
the  derivations  are  Included  In  several  appendices.  Several  additional 
derivations  are  also  Included  In  this  chapter  to  give  the  reader  a feel 
for  what  the  equations  mean  physically. 

Quaslllnear  Equations 

The  equations  of  electrostatic  quaslllnear  theory  governing  the 
evolution  of  the  debris  and  air  distributions  are  derived  from  a simul- 
taneous solution  of  Vlasov's  and  Poisson's  equations  (Ref  11:  27  and  34): 
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The  subscripts  0 and  1 denote  zero  order  and  first  order  quantities 
respectively.  The  equations  are  In  the  CGS-Gausslan  system  of  units. 
These  units  are  used  exclusively  throughout  the  dissertation.  Note  there 
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are  actually  three  Vlasov  equations  In  the  Ion-Ion  Instability  problem. 

One  each  Is  for  the  debris,  air,  and  electrons.  In  the  Vlasov  equations 
for  the  ions,  the  magnetic  field  Is  set  equal  to  zero.  In  the  electron 
Vlasov  equation,  the  field  Is  retained  for  the  cases  where  the  electrons 
are  magnetized.  Otherwise,  the  field  Is  neglected  for  the  electrons  as 
well  because  It  Is  then  Insignificant.  These  latter  cases  are  called  the 
field-free  cases  because  none  of  the  species  Is  Influenced  by  the 
magnetic  field. 

The  Vlasov  equation,  Eqn  (1),  Is  a special  case  of  the  total 
derivative  with  respect  to  time  of  the  distribution  In  Eqn  (3).  This 
distribution  Is  a function  of  the  variables  x,  "v,  and  t.  The  general 
nature  of  the  Vlasov  equation  is  more  apparent  If  the  reader  recalls 
that  velocity  Is  the  time  derivative  of  distance  and  acceleration  Is  the 
time  derivative  of  velocity.  The  acceleration  Is  written  in  terms  of 
the  electric  and  magnetic  fields  causing  the  acceleration.  The  total 
time  derivative  Is  zero  because  the  plasma  Is  colllslonless.  The  distri- 
bution Is  written  Is  two  parts:  The  first  part  Is  the  zero  order, 
spatially  homogeneous  term  which  vprles  very  slowly  with  time.  The 
second  part  Is  the  small,  rapidly  varying  first  order  perturbation  to  the 

I 

zero  order  term.  The  perturbation  Is  caused  by  the  Interaction  of  the 
plasma  particles  with  the  waves,  and  thus,  the  spatial  and  temporal 
variations  of  this  second  term  are  of  the  order  of  the  wavelength  and 
oscillation  period  of  a typical  wave  mode.  Vlasov's  equation  must  be 
solved  simultaneously  with  Maxwell's  equations.  Because  the  waves  are 
electrostatic,  only  Poisson's  equation,  Eqn  (2),  Is  needed. 

Field-Free  Equations.  The  field- free  quaslllnear  equations  are 
derived  following  Drummond  and  Pines  (Ref  12)  and  Bernstein  and 
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Englemann  (Ref  13).  Their  approach  Is  to  spatially  Fourier  transform 
Poisson's  equation  and  the  spatially  averaged,  linearized  Vlasov's  equa- 
tion and  then  solve  the  results  using  a Green's  function.  The  details 
of  this  derivation  are  In  Appendix  A.  An  alternate  derivation  where 
correlation  theory  Is  used  can  be  found  In  Reference  14.  The  resulting 
equations  are 


- o 


= £rl  (io) 


- -i  'n> 


\ (*,*)  = rrr f Jl  .vj,! (12) 

where  >£.  =„^x7?y  wavelength  (13) 

cu(Xj^)  = complex  wave  _ -»  N 


J - vumpieA  wave  , T4  N 

frequency  = ou^(  -t  **  c*>Jl  (14) 

t (X/)  - spectral  wave  x E(-j  S')  E(X,j> 

energy  (ergs)  = v (^Hr  rfr  <15> 

The  Integration  procedure  means  the  Integrals  In  Eqns  (9)  and  (12) 

are  taken  so  the  results  are  analytically  continuous  across  the  oo.  axis 

of  the  complex  wave  frequency  plane.  This  procedure,  equivalent  to 

integrating  along  the  Landau  contour.  Is  covered  in  Reference  15.  This 

reference  also  has  an  excellent  alternate  derivation  of  Eqn  (9). 

The  four  quaslllnear  equations,  Eqns  (9)  to  (12),  form  a set  which 
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must  be  solved  simultaneously  to  evolve  the  Ion  distributions.  Equation 
(9)  Is  the  linear  dispersion  relation  which  determines  the  real  and 

-A 

Imaginary  parts  of  the  waves  corresponding  to  a given  wave  k.  The  real 
part  Is  the  oscillation  frequency,  and  the  imaginary  part  Is  the  growth 
rate.  A positive  Imaginary  part  corresponds  to  a growing  unstable  wave, 
and  a negative  Imaginary  part  corresponds  to  a damped  stable  mode. 
Throughout  the  dissertation,  only  the  fastest  growing  unstable  mode,  If 
any,  at  each  wave  k Is  located  and  used.  Equation  (10)  Is  a diffusion 
equation  In  velocity  space.  This  equation  describes  the  slow  evolution 
of  the  species  distributions  due  to  the  waves.  The  dyadic  diffusion 
coefficient  which  Is  a function  of  time  and  of  the  spatial  variables  Is 
given  by  Eqn  (12).  Note  each  of  the  Ion  species  has  Its  own  separate 
diffusion  equation  and  diffusion  coefficient.  The  diffusion  equation 
for  the  electrons  Is  not  used  because  they  are  adiabatic.  The  only 
coupling  between  the  species  Is  due  to  the  waves  which  are  solutions  of 
the  dispersion  relation.  The  time  rate  of  change  of  the  wave  energy  In 
each  mode  Is  given  by  Eqn  (11). 

The  diffusion  dyad  becomes  very  large  when 

Jx  • ^ ***.iAJ*>)  (16) 

because  > ( -d.  j J')  (17) 

The  particle  velocities  where  Eqn  (16)  holds  define  the  resonant  regions 
of  the  Ion  distributions  because 


-'-Mt  *mim 


ft  a . * 

= wave  phase  velocity  = £:  I svl ^ (18) 

•A 

where  ^ is  the  angle  between  k and  v.  This  resonance  also  dominates 
the  behavior  of  the  dispersion  relation.  The  physical  significance  of 
this  resonance  between  the  wave  phase  velocity  and  the  particle  velocity 
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Is  covered  later  In  the  section  on  Landau  damping/growth.  Note  also 
from  Eqn  (11),  the  wave  mode  energies  change  exponentially  with  time. 
Therefore,  the  diffusion  coefficient  has  a strong  time  depe  idence  via 
the  wave  energies.  These  strong  dependences  of  the  diffusion  coefficient 
on  position  in  velocity  space  and  on  time  make  Eqn  (10)  particularly 
difficult  to  solve. 

The  quasll inear  equations  conserve  matter,  momentum,  and  energy. 

The  proof  of  this  Is  done  by  taking  the  first  three  velocity  moments  of 
Eqn  (10)  as  shown  in  Appendix  B. 

Equations  with  a Magnetic  Field.  Because  the  electrons  are  assumed 
to  be  adiabatic  and  the  Ions  are  assumed  to  be  unmagnetized,  the  only 
quasll Inear  equation  that  changes  when  the  magnetic  field  becomes 
Important  Is  the  dispersion  relation,  Eqn  (9).  This  equation  Is  of  the 
form: 


i - Z = ° 


(19) 


where  is  the  plasma  dielectric  for  species  . The  Ion  dielectrics 
are  unchanged  with  the  addition  of  the  magnetic  field  and  the  electron 
dielectric  becomes  (Ref  16:  158-234) 


*_ 
■*.  2 


( Co  + il.  ) - -H- 


■’■ft  ("  --O.)  ♦ -n.-)£  (*-.). 


(kz  not  near  zero) 
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Co 


or 


_ /x  _ x 


•&/x 


cJ*'-  -tti 


(k2  near  zero)  (20) 


where  to 


= electron 


on  plasma  f 
frequency  = 

(-■O-SX . ? - electron  thermal  speed  [\Tt 

In  the  x,  z direction  =-i/ 

(/  /w* 

= electron  gyro  l *113,1 
frequency  - 


C 


to  + sy\ 


-a. 


/Vi 


llx  sA.  2 Cjt  2 
%(*U)  - plasma  dispersion 


S*\  = 0j  t I 


X - ArA 


function  = 
X x 


“ wk 


cl  i 


2.  - 


X yrl-'*J0LjL 


PLjl  - average  electron  gyro  -J 
(Larmor)  radius  = — 


X'xTjtsr 

syi*-* 


-Jl. 


(21) 

(22) 

(23) 

(24) 

(25) 

(26) 
(27) 


Equation  (20)  is  derived  in  Appendix  C by  averaging  the  electron  motion 

(28) 


in  the  v -v  plane  around  the  magnetic  field  and  by  assuming 
x y 


~ ° 

The  derivation  assumes  the  electrons  are  adiabatic  so  they  remain  In 
their  original  Maxwellian  distribution.  The  well  known  plasma  disper- 
sion function  tabulated  by  Fried  and  Conte  (Ref  17)  results  from  this 
assumption. 

The  magnetized  dispersion  relation  in  the  form  given  by  Eqns  (9) 
and  (20)  Is  valid  when 

C<  / c*l  < < _J1A  (29) 

V/  * 

A x ~ Ion  gyro  radius  = — — > > — > coupling  length  (30) 

'S  -A; 

X < < I (31) 

where  is  the  average  streaming  velocity  of  an  ion  species  and  is 
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an  Ion  gyro  frequency.  The  reasons  for  the  first  two  of  these  conditions 
were  mentioned  in  Chapter  I.  Equation  (31)  insures  the  average  electron 
gyro  radius  is  much  less  than  the  components  of  the  mode  wavelengths  in 
the  x direction  so  the  relatively  simple  form  of  the  dielectric  given  In 
Eqn  (20)  can  be  used. 

Analytical  Form  of  the  Field-Free  Dielectrics.  When  the  ions  are 
in  their  original  Maxwellian  distributions,  an  analytical  form  of  their 
dielectrics  exists.  This  same  form  exists  for  field- free  electrons. 

This  form  results  when  a coordinate  rotation  procedure  is  used  to 
evaluate  the  integral  in  Eqn  (9).  The  dielectrics  are  from  Appendix  D: 


x 

£*■>  ft* 


= \Axcx 


where 


r ***  J 


<*o 

>/x 


7,  (O 


(32) 

(33) 


~ thermal  speed  in  the  x t/x 

wave  direction  = [ (34) 


<**  - \4. 


(35) 


^ a IXI 

The  angle  0 is  that  between  k and  the  x direction.  The  average  velocity 
of  species  is  V^,  and  the  derivative  of  the  plasma  dispersion  function 
defined  in  Eqn  (25)  is 
Anomalous  Transport  Coefficients 

The  anomalous  transport  coefficients  are  obtained  by  taking  velocity 
moments  of  the  spatially  averaged  Vlasov  equation.  (See  Eqn  (A-9).) 

The  standard  hydrodynamic  equations  result  from  this  process  except  that 
anomalous  collision  terms  replace  the  normal  terms  due  to  binary 
collisions.  The  results  from  Appendix  E are  (Ref  18:  1954-1955) 


± f 


V, 
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= - dX  X £(Z,/)  [ 6*  (**>j  X)J  (36) 

—'U 

^■[t^ o'xT0<]  - - x &<  c^j  *£)] 

+ ( u>^  - Jx  • v*  ) I rn  [ £„  ( c*j  X)]^  (37) 

The  right  hand  sides  of  these  equations  are  the  anomalous  collision 
terms.  These  terms  respectively  give  the  momentum  and  energy  density 
transfer  rates  to  the  species  •<  due  to  the  spectrum  of  wave  modes. 
Physical  Weaning  of  the  Quasilinear  Equations 

At  this  time  in  the  dissertation,  a few  discussions  about  the 
physical  meanings  of  the  quasilinear  equations  are  presented.  In  the 
sections  that  follow;  an  analogy  between  velocity  space  and  normal 
diffusion  is  drawn,  the  hydrodynamic  equations  (sometimes  called  the 
fluid  equations)  are  used  to  explain  the  general  behavior  of  waves  in 
the  plasma,  and  a simple  explanation  of  Landau  damping/growth  is  given. 
At  the  end  of  the  last  section,  there  is  a quick  summary  of  the  results 
from  these  discussions. 

Diffusion  Analogy.  An  analogy  can  be  drawn  between  velocity  space 
and  ordinary  diffusion.  In  ordinary  diffusion,  a particle  random  walks 
Its  way  from  one  location  to  another  at  a rate  related  to  the  diffusion 


coefficient: 


0 


where  v Is  the  average  speed  of  the  particle  and  1 is  the  mean  distance 
it  travels  between  collisions  (mean  free  path).  Diffusion  is  a good 
approximation  when  the  mean  free  path  Is  much  less  than  the  characteristic 
size  of  the  system  gradients.  In  velocity  space  diffusion,  a particle 
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random  walks  Its  way  from  one  part  of  the  distribution  to  another  at  a 
rate  related  to  the  diffusion  coefficient: 


% 


K~  <»> 

where  vy  is  the  average  acceleration  and  ly  is  the  velocity  range, av, 
in  which  the  acceleration  takes  place.  Diffusion  is  a good  approximation 
when 


(40) 

where  the  thermal  speed,  cthenn»  is  the  characteristic  size  of  the  grad- 
ients in  velocity  space. 

The  validity  criterion  given  by  Eqn  (40)  can  be  put  in  terms  of  an 
energy  criterion  by  assuming  the  change  in  the  velocity  of  a particle  is 
caused  by  the  absorption  or  emission  of  a quantum  of  wave  energy.  When 
the  wave  energy  is  equated  to  the  change  in  particle  energy, 

W„ave  A Wr^ic/e  = 


- J_ 


/m 


'/>vv  Cf  kfrn. 


(41) 

(42) 

(43) 


where  av  ^ vln1tlftl  and  vinitU1~  ctherm  have  been  assumed.  The 
thermal  energy  of  the  particle  is 

~ (44) 

Therefore,  when  Eqn  (40)  is  multiplied  by  ct^erm  and  Eqns  (43)  and  (44) 
are  used,  the  criterion  becomes 

1VWW<  * < (45) 

Equation  (45)  is  the  weak  turbulence  criterion  of  quasilinear  theory 
implicit  in  the  assumption  that  the  wave  fields  are  first  order 
quantities.  The  term  turbulence  does  not  refer  to  ordinary  spatial 
turbulence.  Rather,  it  refers  to  the  turbulent  random  walk  of  a particle 
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In  velocity  space  caused  by  the  Interaction  between  it  and  many  growing 
unstable  waves. 

Waves  in  an  Unmagnetized  Plasma.  The  oscillation  frequencies  and 
phase  velocities  that  result  from  a solution  of  the  unmagnctized  disper- 
sion relation  are  characteristic  of  the  Interaction  of  the  ions  in  the 
plasma  with  the  electrons.  The  low  frequency  waves  of  interest  charac- 
teristically oscillate  near  the  ion  plasma  frequency  and  propagate  near 
the  Ion  acoustic  speed.  These  properties  can  be  seen  from  a simulta- 
neous solution  of  the  linearized  hydrodynamic  equations  of  continuity, 
momentum  transfer,  and  state  for  each  plasma  species  with  Poisson's 
equation.  The  solution  which  includes  a zero  order  streaming  velocity, 
V^,  is  generalized  from  that  of  Seshadrl  (Ref  7:  98-105).  These  equa- 
tions for  an  initially  homogeneous  plasma  are 


■*Z,]  + v*  ■ 


2 


o 


(46) 


'"do'*'- 


+ fa  ' * XT'' 


E.  =0(47) 


p«, 

'"til 


1 

a* 


(48) 


•V 


Vw  XT- 


(49) 


(50) 


where  a^  = species  <x  sound  speed 
The  constant  Y Is  the  ratio  of  the  specific  heat  capacity  for  constant 
pressure  over  that  for  constant  volume.  When  plane  wave  solutions  of 
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(X'*  - 

the  type  ju"  are  assumed  and  Eqns  (46)  to  (49)  are  combined, 

the  following  dispersion  relation  results: 


i - 7 


x 


- o 


This  is  the  classical  multibeam  dispersion  relation  with  extra  terms 
2 2 

k a^  which  account  for  the  finite  temperatures  of  the  beams. 

In  order  to  look  at  the  characteristic  interaction  between  a single 
ion  species  and  the  electrons,  the  existence  of  only  one  type  of  ions 


V*  9 ° 


are  assumed.  Therefore,  Equation  (51)  becomes 


x 

pA. 


- o 


The  solution  of  this  equation  for  the  low  frequencies  of  interest  is 

y . ycTj,  — ~ i 


'.1  1 
^ JL  j 


The  right  hand  side  of  Eqn  (54)  is  the  general  form  of  the  ion  acoustic 
speed  for  the  ion  species  1.  Thus,  the  wave  propagates  at  the  ion 
acoustic  speed.  If  only  one  ion  species  exists,  it  is  singly  ionized  so 


x 

'.U  ^ 


= / 


and  Te>>  T^  the  often  quoted  expression  for  this  speed  results: 

■ i 


The  speed  is  a function  of  the  electron  temperature  and  the  ion  mass  and 
is  much  slower  than  the  electron  sound  speed  (Eqn  (50)  with*  = e).  The 
phase  velocity  is  slow  because  the  light,  rapidly  moving  electrons 


24 


T — 

l 


AFIT/DS/PH/78-1 


generate  electric  fields  which  can  only  slowly  transfer  the  electron 

momentum  to  the  heavy  ions.  (A  more  rigorous  explanation  is  in  Reference 

7:  102-103.)  The  condition  T ^>T.  assures  the  acoustic  waves  are  not 

e 1 

strongly  Landau  damped,  a process  which  is  described  later  in  this 
chapter. 

The  oscillation  frequency  is  obtained  directly  from  Eqn  (53)  by 
assuming  co <<  60pg  and  a^  = 0 (T^  <<Tg).  The  result  is 


■x  _ 


pji. 

I + ^ 


■z 


^~D  JS. 


where  2 n = electron  Debye  length 


_o  -C 


Therefore,  when  the  wavelength  is  less  than  the  electron  Debye  length, 
the  solution,  Eqn  (57),  is 

co  2=  (60) 

When  the  reverse  holds,  the  oscillation  frequency  Is  somewhat  below  the 
ion  plasma  frequency.  These  two  cases  respectively  correspond  to 
stationary  (Vp^  = 0)  electrostatic  ion  oscillations  and  to  the  ion 
acoustic  wave  with  a phase  velocity  given  by  Eqn  (56).  To  summarize  the 
above  results  for  waves  in  an  unmagnetized  plasma,  the  characteristic 
wave  phase  velocity  is  near  the  ion  acoustic  speed,  and  the  character- 
istic oscillation  frequency  is  near  or  just  below  the  ion  plasma 
frequency. 

Waves  in  a Magnetized  Plasma.  The  general  nature  of  waves  in  a 
magnetized  plasma  can  be  seen  from  fully  electromagnetic  magnetoionic 
theory  (Ref  16:  9-12).  Kagnetoionic  theory  results  from  the  simulta- 
neous solution  of  the  hydrodynamic  equation  of  momentum  transfer,  Eqn 
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(47),  with  Maxwell's  equations.  The  theory  assumes  a cold  plasma  which 
means  that  the  pressure  term  Is  set  equal  to  zero.  Also,  no  zero  order 
drift  velocities  exist.  In  the  derivation  below,  the  waves  are  assumed 
to  propagate  In  the  x direction  perpendicular  to  a uniform  magnetic 
field  which  Is  In  the  z direction. 

The  type  of  wave  which  Is  of  Interest  Is  the  transverse  magnetic 
(TM)  wave.  This  wave  has  its  magnetic  field  parallel  to  the  ambient 
magnetic  field.  Its  electric  field  lies  In  the  x-y  plane  as  shown  In 
Fig.  2.  If  the  electric  field  Is  very  nearly  In  the  propagation 


B B 
o wave 


/ / wave 
x I (in  x-y  plane) 

x " 

Fig.  2.  Transverse  Magnetic  Wave. 

direction,  the  TM  mode  becomes  nearly  electrostatic.  The  dispersion 
relation  for  this  type  of  wave  Is  (Ref  7:  328) 


~ _ f e, ' 

C-A~  ~ / ce.  + off.-o. 


top* 


z., 

C*Jp^ 


where  (,  = I - ^ ~ tjZZ 


£%  ' slX) 


Again,  the  existence  of  only  one  species  of  Ions  has  been  assumed. 
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Equation  (61)  can  be  rewritten  as: 

« ,/'i 
' [ (c*  - - S3%>] 


where  co  = normalized  oscillation 

frequency 

Z3\  - normalized  cut  off  _ 
x freouencles  = x . 


(67) 

/4 


co i = normalized  cut  off  _ f ^ 

x frequencies  = 4-  + 'J  {66) 

wuH  = normalized  upper  hybrid  f V/a 

resonant  frequency  - I I + — j (67) 

coLH  = normalized  lower  hybrid  r ^ ^ j J4 

resonant  frequency  j (68) 

Note  the  frequencies  In  Eqns  (65)  to  (68)  are  normalized  to  the  electron 

plasma  frequency.  The  unnormalized  form  of  the  lower  hybrid  frequency  Is 

coLH  - yi  +.  *>?*/. nl  ^69) 

The  frequency  range  for  Ion  wave  modes  satisfies: 

cZx  <<  <7°) 


so  Eqn  (64)  becomes 


This  equation  is  valid  for  the  frequency  range: 

O 4 co  < &Oi_H 

The  limit  of  Eqn  (71)  as  the  frequency  goes  to  zero  is 


* '/yJ~aX 


where  v#lf  - Alfven  wave  speed  = W 

The  ambient  field  is  Bq.  Therefore,  the  characteristic  velocity  of  the 
TM  mode  Is  the  Alfven  speed.  The  mode  propagates  approximately  at  this 
speed  except  when  the  frequency  Is  almost  exactly  at  the  lower  hybrid 
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frequency  (Ref  7:  330). 

The  transverse  magnetic  mode  is  largely  electrostatic  when  (Ref  7: 


400-402) 


l^irl  »»  It.IJiJ 


The  electrostatic  version  of  Eqn  (61), 


^ i ' I cJ*' — cox~  -fiJJ  ” ” 
results  when  Eqn  (75)  Is  assumed.  The  solution  of  this  equation  in  the 
frequency  range  given  by  Eqn  (29)  Is 
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ur-  -n. 


= o 


6o  - 
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1/1  ^/Cn.2 

From  Eqns  (72)  and  (77),  the  characteristic  frequency  of  oscillation  of 
the  TM  mode  Is  close  to  but  somewhat  below  the  lower  hybrid  frequency. 
The  electrostatic  criterion,  Eqn  (75),  Is  approximately  equivalent 


to  the  criterion: 


OJp*. 

SI.  3 


This  equation  Is  used  to  test  the  validity  of  the  electrostatic  approx- 
imation for  the  waves  in  a magnetized  plasma. 

The  lower  hybrid  frequency  has  an  interesting  physical  significance. 
The  high  density  limit  of  Eqn  (69),_n.e/co  pe<< 


1/Sl-  -ft-j 


In  this  limit,  the  lower  hybrid  frequency  Is  the  geometric  mean  of  the 
electron  and  Ion  gyro  frequencies.  This  same  result  can  be  obtained  by 
setting  the  net  electron  polarization  drift  velocity  across  the  magnetic 
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field  equal  to  the  Ion  velocity.  The  electron  velocity  Is  (Ref  7:  201) 
d . p 

sr  - ^£l_  - rr,^c  C eo  < < -O  (80) 

^ 60_a*  5.^* 

where  harmonic  time  dependence  for  the  electric  field  Is  assumed.  The 
Ion  velocity  from  Newton's  second  law, 

^ (81) 


For  simplicity,  the  ions  are  assumed  to  be  singly  charged.  When  the 
velocities  In  Eqns  (80)  and  (82)  are  equated,  the  lower  hybrid  frequency 
in  Eqn  (79)  results.  At  this  frequency,  the  electrons  and  ions  drift 
together  so  the  plasma  has  no  net  current. 

Unmagnetlzed  vs  Magnetized  Haves.  The  usual  method  of  determining 
whether  the  plasma  is  magnetized  or  not  is  by  its  "beta"  which  is  the 
ratio  of  the  particle  pressure  to  the  magnetic  field  pressure.  Because 
only  the  electrons  may  be  magnetized,  the  relevant  beta  is  that  for  the 


electrons: 


. _ thermal  pressure  = ^ — (B3, 

magnetic  pressure  * B*/  t'Tt'  V 

where  <<  I magnetized  (84) 

or  > > l unmagnetlzed  (85] 

The  beta  also  relates  the  relative  importance  of  the  ion  acoustic  speed 

and  the  Alfven  wave  speed: 


3 ^ 
C 4. 


where  the  speeds  are  given  by  Eqns  (56)  and  (74).  The  range  of  validity 
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of  the  magnetized  equation,  however.  Is  determined  by  the  parameter 
given  by  Eqn  (26).  The  relationship  between  beta  and  X Is 


X * X ( l o' 


where  is  the  relative  velocity  In  cm/sec  and  AMU  Is  the  debris  atomic 
weight  In  atomic  mass  units.  In  the  derivation  of  Eqn  (87),  the  following 


were  assumed: 
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Equation  (88)  states  the  phase  velocity  Is  of  the  order  of  the  relative 

velocity.  The  reason  for  this  will  become  apparent  In  the  next  section. 

Equation  (89)  Is  true  because  the  oscillation  frequency  is  near  the  ton 

plasma  frequency  or  the  lower  hybrid  frequency.  Because  the  relative 

8 

velocities  of  interest  are  near  10  cm/sec,  the  magnetized  equation 
validity  criterion,  Eqn  (31),  approximately  corresponds  to  Eqn  (84). 
Therefore,  the  choice  of  the  magnetized  versus  the  unmagnetized  approach 
Is  consistent  with  the  use  of  the  magnetized  equations  when  Eqn  (31) 
holds  and  the  use  of  the  unmagnetized  set  when  It  does  not. 

Landau  Damping/Growth.  The  growth  rates  of  the  modes  in  the  un- 
stable wave  spectrum  are  mostly  determined  by  a resonant  Interaction 
between  planes  of  constant  wave  phase  and  an  unstable  part  of  the  Ion 
distributions.  The  unstable  part  for  plasmas  where  the  two  Ion  beams 
have  a moderately  low  relative  velocity  Is  shown  In  Fig.  3.  The  figure 
Is  a top  view  of  the  two-dimensional  debris  and  air  distributions 
plotted  as  contours  In  the  velocity  plane.  (The  reference  frame  In  this 
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Fig.  3.  Resonance  Between  the  Unstable  Region  and  an  Unstable  Kave. 
figure  Is  obtained  using  a procedure  described  in  the  next  chapter.)  The 
distance  between  the  peaks  of  the  debris  and  air  "mountains"  is  the 
relative  velocity  between  the  beams  (the  expansion  velocity).  A typical 
wave  mode  propagating  In  the  vx-v2  plane  Is  shown  with  Its  phase  planes 
moving  In  step  with  the  unstable  part.  This  resonance  between  the 
particles  In  the  unstable  part  of  the  distributions  and  the  waves  Is  the 
resonance  given  by  Eqn  (18).  Therefore,  the  unstable  part  Is  the 
resonance  region  of  the  Ion  distributions.  Because  the  region  Is  near 
the  vx  axis,  the  angles  ^ and  9 are  approximately  equal.  Also,  the 
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resonant  particles  have  velocities  which  are  approximately  equal  to  the 
approximately  constant  relative  velocity.  Thus, 

~ Vn  tonj+«nt  (90) 

C-rQ.  & * 

The  resonant  region  Is  unstable  because  It  has  a positive  slope 
meaning  that  more  plasma  Ions  are  at  higher  velocities  than  are  at  lower 
velocities.  The  preferred  equilibrium  state  of  the  plasma  Is  a single 
combined  Ion  distribution  with  more  Ions  at  lower  velocities  than  at 
higher  ones.  Therefore,  as  the  plasma  relaxes  toward  Its  equilibrium 
state,  the  excess  particle  kinetic  energy  must  be  given  to  the  growing 
waves.  This  Landau  growth  effect  Is  the  opposite  of  the  usual  Landau 
damping  process  which  occurs  for  stable  plasmas  (Ref  19:  383-392). 

An  Interesting  two-dimensional  projection  effect  occurs  with  the 
unstable  waves.  From  Eqn  (90),  only  the  phase  velocity  over  cos©  has  to 
match  the  relative  velocity.  If  the  phase  velocity  Is  much  smaller  than 
the  relative  velocity,  the  resonance  can  still  occur  If  the  angle  ©Is 
close  to  90  degrees.  This  projection  effect,  shown  In  Fig.  3,  allows  a 
wave  propagating  at  the  Ion  acoustic  speed  or  the  Alfven  speed  to 
resonate  with  a pair  of  Ion  distributions  having  a much  larger  relative 
velocity.  As  a matter  of  fact,  the  terms  "low  relative  velocity  plasma" 
and  "high  relative  velocity  plasma"  which  are  used  throughout  the  dis- 
sertation refer  to  this  effect.  In  a low  velocity  plasma,  the  relative 
velocity  Is  about  the  same  as  the  characteristic  wave  phase  velocity: 
either  the  Ion  acoustic  speed  or  the  Alfven  speed.  The  unstable  waves 
for  this  type  of  plasma  propagate  nearly  In  the  streaming  direction.  In 
a high  velocity  plasma,  the  relative  velocity  is  much  larger  than  the 
characteristic  wave  velocity  so  the  waves  propagate  at  large  angles  to 
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the  streaming  direction. 


The  exact  resonant  regions  for  large  relative  velocity  plasmas  lie 
somewhat  off  the  v axis  as  shown  In  Fig.  4,  a close  up  view  of  the  air 


PHASE  PLANE 


RESONANT  REGIONS 


Efield 

Fig.  4.  Exact  Resonant  Regions. 

distribution.  The  resonance  between  the  typical  wave  mode  and  the  reso- 
nant region  below  the  axis  Is  shown.  The  equivalent  wave  mode  propa- 
gating in  the  -B  direction  resonates  with  the  region  above  the  axis. 

(The  symmetries  of  wave  modes  propagating  in  various  directions  is 
discussed  In  Appendix  F.  See  Eqns  (F-39)  to  (F-42).)  Because  the  wave 
modes  are  electrostatic,  their  oscillating  electric  fields  can  be  expected 
to  diffuse  the. ions  primarily  in  the  wave  propagation  directions.  Hence, 
the  resonant  regions  lie  approximately  at  the  same  angle  to  the  vx  axis 
as  the  fastest  growing  mode  angle. 

The  unstable  waves  can  only  appear  if  the  electrons  do  not  freely 
move  around  to  cancel  out  the  electric  fields  generated  by  the  ion 
motions.  The  electrons  can  be  hindered  by  a high  temperature  or  by  the 
applied  magnetic  field.  Hot  electrons  have  a relatively  large  Debye 
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length  (See  Eqn  (59).)  so  their  shielding  effect  on  the  ion  generated 
fields  is  reduced.  Therefore,  for  the  unmagnetized  cases,  unstable  wave 
modes  can  be  expected  for  plasmas  where  the  electron  Debye  length  is 
comparable  to  the  mode  wavelengths.  For  the  magnetized  cases,  unstable 
modes  propagating  approximately  perpendicular  to  the  magnetic  field  can 
exist  for  a much  colder  electron  distribution.  This  occurs  because  the 
field  inhibits  the  electron  motion  perpendicular  to  it.  The  angular 
range  of  these  modes  around  the  streaming  direction  would  be  relatively 
small  because  the  cold  electrons  can  freely  move  along  the  magnetic  field. 
The  allowed  free  motion  of  the  electrons  along  the  field  also  explains 
why  electromagnetic  effects  tend  to  reduce  the  growth  rates  of  waves  in 
a magnetized  plasma.  The  wave  magnetic  field  of  an  electromagnetic  wave 
bends  the  ambient  magnetic  field  allowing  the  electrons  to  move  more 
easily  along  it.  (See  also  Ref  20.) 

As  a summary,  several  general  conclusions  can  be  drawn  from  the 
discussions  in  the  above  sections.  They  are  as  follows: 

(1)  Equation  (10)  is  a diffusion  equation  in  velocity  space. 

(2)  The  wave  mode  oscillation  frequencies  are  determined  by  the 
plasma  parameters  such  as  density  and  magnetic  field  strength. 

(3)  The  mode  wavelengths  for  a given  plasma  are  determined  by  the 
relative  velocity  between  the  debris  and  the  air. 

(4)  The  growth  rates  of  the  wave  modes  are  determined  by  the  shape 
of  the  ion  distributions  and  the  relative  mobility  of  the  electrons. 
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III . The  Numerical  Solution  Procedure 

t 

In  this  chapter,  the  general  solution  procedure  for  the  quaslllnear 
equations  is  described.  First,  the  reference  frame  chosen  for  the 
solution  is  presented.  Next,  the  procedure  used  to  solve  the  equations 
is  described.  Then,  the  recasting  of  the  equations  into  a form  suitable 
for  numerical  solution — by  no  means  a trivial  process--is  covered.  And 
finally,  the  computer  program  is  described.  As  before,  many  of  the 
details  are  relegated  to  the  appendices. 

Before  these  discussions  are  started,  an  important  notational 
distinction  must  be  made:  From  this  point  on  in  the  dissertation,  the 
term  "distribution"  can  either  refer  to  the  individual  species  velocity 
distributions  as  before,  or  it  can  refer  to  the  type  of  plasma  as  a 
whole.  The  appropriate  meaning  will  always  be  clear  from  the  context. 
Usually,  the  word  "distribution"  (singular)  will  refer  to  the  type  of 
plasma,  and  the  word  "distributions"  (plural)  will  refer  to  the  plasma 
species. 

Plasma  Distribution 

A typical  assumed  initital  plasma  distribution  in  an  earth  station- 
ary reference  frame  is  shown  in  Fig.  5.  The  figure  shows  the  nuclear 
debris  ions  and  the  air  ions  superimposed  on  the  hot  electron  background. 
The  top  half  is  a side  view  in  the  f*-vx  plane,  and  the  bottom  half  is 
the  top  view  showing  the  two  dimensionality  of  the  species  distributions. 
Again,  the  electrons  come  from  both  ion  species  and  are  premixed  by  fast 
processes  that  would  be  essentially  complete  before  the  ion-ion  evolution 
begins.  All  three  distributions  are  Maxwellian.  The  average  velocity 
of  the  air  is  zero,  and  the  relative  velocity  is  the  speed  the  debris 
moves  through  the  air.  The  plasma  distribution  is  spatially  homogeneous 
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Fig.  5.  Typical  Plasma  Distribution  In  an  Earth  Stationary 
Reference  Frame. 
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because  It  represents  the  plasma  In  the  small  volume  In  the  debris-air 
mixing  region  shown  In  Fig.  1. 

Figure  6 shows  the  same  distribution  In  the  reference  frame  used  In 
the  computer  calculations.  This  figure  has  its  vx  axis  reversed  from  the 
previous  figure.  The  average  electron  velocity  and  net  current  are  set 
equal  to  zero.  The  debris  plasma  frequency  Is  assumed  to  be  greater 
than  the  air  plasma  frequency  so  the  fastest  growing  wave  modes  are 
directed  toward  the  right.  This  choice. is  for  convenience  only  because 
the  roles  of  the  "debris"  and  "air"  beams  can  be  easily  interchanged. 

The  phase  velocity  of  a typical  wave  mode  propagating  at  an  angle  & to 
the  streaming  direction  Is  shown.  The  magnetic  field  is  uniform  and  is 
In  the  z direction. 

This  reference  frame  is  used  In  the  solution  of  the  quasllinear 
equations  because  a known,  time  independent,  zero  order  solution  must 
exist  for  the  linearization  procedures  described  in  the  previous  chapter 
to  be  valid.  In  this  frame,  the  electrons  have  no  drift  velocity 

U — x 

relative  to  the  magnetic  field  so  there  Is  no  zero  order  V X B motion. 
Because  the  Ions  are  unmagnetized  as  explained  earlier,  their  zero  order 
drift  Is  purely  In  the  x direction,  and  no  V X B motion  exists  for  them 
either.  The  Ion  drift  velocities  are  arranged  to  assure  zero  order 
charge  neutrality,  to  set  their  difference  equal  to  the  desired  relative 
velocity,  and  to  set  the  net  current  equal  to  zero.  The  zero  current 
assures  that  there  Is  no  zero  order  J X B motion  of  the  plasma. 

An  Important  point  to  note  at  this  time  Is  that  the  plasma  para- 
meters stated  In  the  dissertation  are  given  In  this  reference  frame. 
Again,  this  Is  the  frame  shown  in  Fig.  6.  Examples  of  such  parameters 
are  the  species  streaming  (average)  velocities  and  the  wave  phase 
velocities. 
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*7 


Calculatlve  Procedure 

The  basic  calculatlve  procedure,  which  was  outlined  in  Chapter  I,  is 
to  trace  the  evolution  in  time  of  the  ion  distributions  as  they  interact 
with  a spectrum  of  growing  waves.  This  evolution  Is  the  solution  of  the 
quasilinear  equations.  The  distributions  are  followed  in  time  until  the 
anomalous  transport  coefficients  begin  declining  indicating  the  approach 
to  saturation.  The  calculation  becomes  prohibitively  expensive  when 
carried  beyond  this  point.  In  an  evolution,  the  energy  in  each  unstable 
mode  in  the  wave  spectrum  is  followed  from  its  original  thermal  back- 
ground value  to  its  final  value  when  the  calculation  terminates.  The 
calculation  proceeds  by  alternately  computing  the  waves  and  evolving  the 
ion  distributions.  This  procedure  makes  use  of  the  quasilinear  assumption 
that  the  wave  frequencies  and  growth  rates  change  relatively  slowly. 

The  evolution  calculation  Is  started  by  finding  the  unstable  wave 
modes.  If  any,  from  the  initial  Maxwellian  species  distributions.  These 
solutions  of  the  linear  problem  use  the  analytical  forms  of  the  dispersion 
relation  because  the  results  can  be  generated  more  quickly  and  accurately. 
The  linear  problem  must  be  solved  many  times  because  the  program  must 
locate  the  wavelength  and  propagation  direction  range  of  the  unstable 
modes  and  then  compute  those  modes  which  are  followed  during  the  evolu- 
tion of  the  ions. 

The  ion  evolution  begins  when  the  wave  spectrum  Is  allowed  to  act 
on  the  Ion  distributions  causing  them  to  change  In  time.  This  evolution 
Is  that  described  by  the  quasilinear  diffusion  equation.  The  evolution 
is  continued  until  enough  changes  In  the  Ion  distributions  have 
accumulated  so  that  the  wave  frequency  and  growth  rate  of  each  mode  In 
the  spectrum  must  be  recomputed.  The  growth  rate,  which  changes  faster 
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than  the  frequency,  is  a function  of  the  ion  distribution  slopes.  Thus, 
the  program  uses  the  slope  of  the  faster  changing,  less  dense  air  distri- 
bution to  determine  when  to  temporarily  stop  the  evolution  of  the  ions  In 
order  to  recompute  the  waves.  (The  mode  update  procedure  will  be 
explained  in  more  detail  in  the  description  of  the  computer  program.) 

Because  the  recomputation  of  the  waves  uses  evolved  ion  distributions 
which  are  no  longer  Maxwellian,  a numerical  form  of  the  dispersion 
relation  integrals  must  be  used.  The  computation  of  these  integrals  is 
relatively  slow  so  the  modes  have  to  be  updated  as  Infrequently  as 
possible.  The  updated  modes  have  smaller  growth  rates  because  the 
general  effect  of  the  changes  In  the  ion  distributions  is  to  stabilize. 

The  anomalous  transport  coefficients  are  computed  just  after  the  modes 
are  updated  because  the  integrals  for  these  coefficients  need  the  current 
values  of  the  frequencies  and  growth  rates. 

The  calculation  continues  by  alternately  evolving  the  ion  distri- 
butions and  updating  the  wave  modes  until  it  is  finished.  Each  cycle  of 
wave  mode  update  and  ion  evolution  is  an  "evolution  phase."  The 
anomalous  transport  coefficients  which  are  computed  after  each  mode  up- 
date initially  grow  in  time  as  the  wave  energy  in  the  spectrum  increases. 
When  the  evolution  begins  to  saturate,  the  coefficients  peak  and  then 
begin  to  decline  because  the  mode  growth  rates  and  the  imaginary  parts 
of  the  ion  plasma  dielectrics  are  decreasing  faster  than  the  wave  energy 
Is  increasing.  (See  Eqns  (36)  and  (37).) 

Solution  Program 

The  computer  program  written  to  solve  the  four  quasi! inear  equations 
(Eqns  (9)  to  (12))  is  based  on  a set  of  routines  In  the  GEARB  package 
from  the  Lawrence  Livermore  Laboratory.  The  GEARB  package  solves  a set 
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of  simulaneous  ordinary  differential  equations  resulting  from  a finite 
differencing  of  the  diffusion  equation,  Eqn  (10),  by  the  method  of  lines 
(Refs  21,  22,  and  23).  The  coefficients  for  this  set  of  equations 
include  the  diffusion  coefficient,  Eqn  (12),  which  is  recast  in  a suitable 
numerical  form.  The  other  two  equations,  Eqns  (9)  and  (11),  are  used  to 
calculate  the  wave  mode  frequencies,  growth  rates,  and  spectral  energies 
which  are  needed  in  the  diffusion  coefficient  integrals. 

The  computer  solution  takes  place  on  three  separate  levels:  (1)  the 
complex  frequency  plane  on  which  the  solutions  to  the  dispersion  relation 
are  found,  (2)  the  kx-k2  plane  on  which  the  unstable  waves  are  located 
and  the  diffusion  coefficient  integrals  are  done,  and  (3)  the  vx-v2  plane 
on  which  the  ion  distributions  are  evolved  and  the  dispersion  relation 
integrals  are  done.  Actually,  there  are  three  separate  velocity  planes, 
one  each  for  the  debris,  air,  and  electrons. 

Recasting  the  Quasi! inear  Equations.  The  details  of  the  recasting 
of  the  quasilinear  equations  are  in  Appendix  F.  The  general  recasting 
procedure  begins  with  a transformation  of  the  entire  velocity  plane  onto 
a square.  The  spatial  variables  of  the  diffusion  equation  in  the  trans- 
formed coordinates  are  then  finite  differenced  by  the  method  of  lines. 
Also,  the  numerical  forms  of  the  diffusion  coefficient  integrals  and  of 
the  dispersion  relation  integrals  are  constructed.  The  dispersion 
relation  integrals  are  analytically  continued  so  the  dispersion  relation 
results  are  valid  for  damped  [o>^  < 0)  modes.  The  analytical  continuation 
of  the  diffusion  coefficient  integrals  is  not  done  because  the  evolutions 
are  always  finished  when  most  of  the  modes  still  have  positive  growth 
rates.  The  programming  of  the  analytical  continuation  of  these  integrals 
is  a very  difficult  task. 
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Computer  Program.  The  quaslllnear  equations  are  solved  by  a set  of 
computer  routines  controlled  by  a master  program,  IONION.  This  program 
sets  up  the  problem,  controls  the  solution  of  the  linear  problem, 
controls  the  evolution  of  the  ion  distributions  in  the  manner  described 
earlier,  and  outputs  the  results. 

Several  sets  of  routines  were  obtained  to  perform  specific  tasks. 

The  first  of  these,  the  GEARB  package  mentioned  above,  solves  the  coupled 
ordinary  differential  equations.  The  second  set  of  routines,  CONTOUR,  is 
a complex  root  finder  (Ref  24).  This  set  Is  used  to  find  the  complex 
frequency  solutions  of  the  dispersion  relation  for  the  initial  plasma 
distribution.  The  next  routine,  CROOT,  Is  used  throughout  the  evolution 
to  Improve  dispersion  relation  solutions.  The  fourth  set  of  routines; 
SINGRT,  SINTRP,  and  XYZSPLN;  does  differentiations  and  Integrations.  The 
fifth  set,  DISP,  Is  used  to  numerically  generate  the  values  of  the 
plasma  dispersion  function  and  its  derivative.  The  routines  CROOT 
through  DISP  are  available  from  the  AF1JL  Math  Library.  The  sixth  group 
of  routines;  SAVEIT,  BOSS,  and  SNAP;  was  obtained  from  Harry  Murphy  at 
AFWL/DYVM.  These  are  utility  routines  to  write  restart  files,  time 
functions,  and  facilitate  debugging.  A final  set  of  routines,  available 
on  the  AFWL  CDC  computer  system  library,  generates  plots. 

Most  of  the  routines  listed  above  were  modified  before  Incorporation 
Into  the  solution  program.  A more  complete  explanation  of  these  routines 
and  of  those  written  for  the  dissertation  Is  in  Appendix  G.  The  appendix 
also  contains  a schematic  diagram  to  show  how  the  solution  program  Is 
organized.  A complete  listing  is  available  from  the  dissertation  author 
at  AFVL/DYCE,  Klrtland  AFB,  New  Mexico  87117.  Most  of  the  routines  are 
written  In  Fortran.  The  total  program  length  Is  about  25,000  decimal 
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words . 


An  Interesting  feature  of  the  solution  program  Is  that  It  Is  also 
set  up  to  evolve  the  Ions  and  electrons  for  magnetic- field- free, 
electron-ion  Instabilities.  This  feature  has  never  been  tested,  however. 
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IV.  The  Computational  Results 

In  this  chapter,  the  results  from  the  computer  runs  are  described. 

In  the  first  few  sections  that  follow,  the  tests  of  the  linear  solution 
and  the  test  of  the  full  evolution  are  discussed.  Several  test  cases 
were  done  and  were  compared  with  analytical  solutions  to  insure  the 
program  was  working  properly.  Only  the  major  tests  of  the  program  as  a 
whole  are  described.  Each  subroutine  was  also  tested  before  it  was 
incorporated  into  the  main  program.  Next  in  this  chapter,  a few  of  the 
results  of  debris-air  coupling  problem  simulations  at  200  km  and  600  km 
are  presented.  Only  a few  highlights  are  given  to  avoid  inundating  the 
reader  with  data.  A general  discussion  of  these  latter  results  is  in 
Chapter  V. 

Tests  of  the  Linear  Dispersion  Relation  Solution 

The  computer  results  for  the  linear  problem  (the  dispersion  relation 
solutions  for  the  initial  Maxwellian  distributions)  were  tested  using 
three  different  types  of  plasmas  for  which  analytical  solutions  are 
relatively  easy  to  obtain.  The  analytical  and  computer  linear  solutions 
for  each  of  these  types  are  compared  in  the  next  three  sections. 

Equal  and  Opposite  Cold  Ion  Beams.  The  first  test  case  is  for 
equally  cold,  equally  dense  debris  and  air  distributions.  The  term 
"cold"  means  the  thermal  widths  of  the  ion  distributions  are  much  less 
than  the  relative  velocity  between  the  beams.  The  ion  atomic  weights 
are  also  assumed  to  be  the  same  so  in  the  computer  reference  frame  the 
ion  bears  have  equal  and  opposite  velocities.  The  unmagnetized  electrons 
are  assumed  to  be  hot  so  unstable  wave  modes  can  exist.  The  relative 
velocity  is  chosen  to  be  less  than  the  ion  acoustic  velocity  so  the 
problem  is  essentially  one  dimensional.  Because  the  problem  is  one 
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dimensional,  unstable  wave  modes  which  propagate  In  the  streaming  direc- 
tion can  exist.  The  analytical  solution,  therefore,  Is  done  In  one 
dimension  and  Is  compared  with  the  computer  generated  modes  where  kx  «=  k 
and  kz  = 0. 

The  Ions  are  cold  so  the  wave  phase  velocities  are  well  away  from 
the  bulk  of  their  distributions: 

— - V*  | > C*  * - ».  * (91) 

Also,  because  the  electrons  are  hot, 

|3r|  <<  (92) 


Therefore,  the  problem  Is  hydrodynamic,  and  the  dispersion  relation  In 
Eqn  (51)  applies.  In  hydrodynamic  problems,  the  effect  of  the  distribu- 
tion shapes  can  be  ignored.  This  case  Is  very  similar  to  the  equal  and 
opposite  electron  beam  problem  frequently  found  In  texts.  (For  example, 
see  Ref  19:  222-227.) 


Equation  (51)  after  Eqns  (91)  and  (92)  are  used  becomes 


1 + 

•x  . 

— ^JOjC 

f ' + * - - 1 

(93) 

Jfcl 

where 

(94) 

v = 

ll 

0 

1 

(95) 

--  c*  x.  - 1 . y.  - k,  * 3 (96) 


Note  the  value  for  Ye  from  Ref  25:  189  Is  the  value  traditionally  used. 
Equation  (93)  also  results  when  a power  series  expansion  procedure  Is 
applied  directly  to  the  non  hydrodynamic  dispersion  relation,  Eqns  (19) 
and  (32)  (Ref  19:  381-392). 

The  solution  of  Eqn  (93)  Is  (Ref  26:  129-131) 
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= \ax,a^)  [u^JE^  * ux+^&)J  (97) 


where  r 


(98) 


D.* 


Unstable  wave  modes  exist  when  Imaginary  solutions  for  the  In  Eqn  (97) 
exist.  Hence  for  Instability,  the  minus  sign  must  be  used,  and  the  wave 
k must  be  In  the  following  range: 

I I 


° 4 4 x ( v* 


(99) 


where  v is  given  by  Eqn  (56).  Because  k Is  real, 

dC 

>wu*f 

v < ^ 


(loo) 


for  Instability.  Equation  (100)  is  the  usual  one-dimensional  instability 
criterion.  The  solution  of  Eqn  (97)  for  a wave  k which  satisfies  Eqn 
(99)  is  the  imaginary  part  of  the  complex  frequency  for  the  growing  mode 
at  that  k.  The  real  part  is  always  zero. 

The  test  calculations  are  done  using  the  plasma  shown  in  Table  H-I 
in  Appendix  H.  Tabulated  in  this  appendix  are  all  of  the  types  of 
plasmas  used  for  the  dissertation  calculations.  The  beam  densities  are 
not  made  exactly  equal  because  that  would  cause  the  unstable  roots  to 
lie  on  the  imaginary  axis  of  the  complex  frequency  plane.  The  plasma 
dispersion  function  varies  many  orders  of  magnitude  across  that  axis  so 
the  program  gives  bad  results  due  to  numerical  errors.  The  values  of 
the  ion  densities  are  chosen  to  approximate  those  found  at  200  km  above 
the  earth  if  the  air  there  were  completely  singly  ionized  (Ref  27). 

For  the  distribution  in  Table  H-I,  the  unstable  wave  k range  from 
Eqn  (99)  is 
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o < A < <i.i  «**-'  Uoi) 

and  the  complex  frequency  solution  of  Eqn  (97)  for  k * 6.0  cm-*  is 

Caj  = o.o  -#■  l.o  ( to'*)  Sec-'  (102) 

The  computer  solutions  agree  well  with  the  results  In  Eqns  (101)  and 
(102).  The  unstable  wave  k range  Is  from  2 cm-*  to  9 cm-*,  and  the 
complex  frequency  for  k s 6 cm"*  is  ^>  = 8.6(10^)  + 1 3.0(10^)  sec-*.  The 
small  real  part  In  the  wave  frequency  is  due  to  the  slight  density 
inequality.  Unstable  modes  are  not  found  for  wave  k's  less  than  2 cm-* 
due  to  the  following  thermal  effect:  The  phase  velocity  of  each  unstable 
mode  must  match  the  velocity  of  the  distribution  resonant  region  shown 
in  Fig.  3.  Yet,  the  oscillation  frequency  must  remain  near  the  character- 
istic one  which,  in  this  case,  is  determined  by  the  ion  plasma 
frequencies.  From  the  definition  of  the  phase  velocity  in  Eqn  (18),  a 


small  wave  k implies  too  small  of  oscillation  frequency  so  the  plasma 
waves  cannot  resonate.  This  effect  is  absent  in  the  cold  plasma  hydro- 
dynamic  equations  because  the  distributions  are  assumed  to  be  infinitely 
narrow. 

An  Interesting  note  is  that  the  computer  solutions  found  numerous 
unstable  wave  modes  which  propagate  at  large  angles  to  the  streaming 
direction.  Unstable  modes  existed  at  angles  from  zero  degrees  to  75 
degrees.  This  result  is  not  predicted  by  one-dimensional  theory. 

Bump-on- the-Tail  Distribution.  The  next  test  case  is  the  bump-on- 
the-tail  distribution  which  is  a weak,  cold  air  beam  passing  through  a 
warm  debris  plasma.  Like  the  equal  beam  distribution,  this  one  also  has 
a commonly  done  electron  analogy  (Ref  19:  458-463).  Again,  the  atomic 
weights  of  the  air  and  the  debris  are  assumed  to  be  the  same.  The  debris 
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beam  Is  assumed  to  be  warm  and  much  denser  than  the  air  beam.  The  air 
beam  is  assumed  to  be  cool  and  to  have  a mean  velocity  much  greater  than 
the  thermal  velocity  of  the  debris.  Again,  the  electrons  are  assumed  to 
be  hot  and  unmagnetized.  A side  view,  velocity  space  plot  of  the  sum  of 
the  debris  and  air  distributions  In  the  computer  reference  frame  Is 
shown  In  Fig.  7.  The  initial  sum  (the  solid  line  in  the  figure)  Is  that 


Fig.  7.  Bump-on-the-Tall  Distribution  Before  and  After  Evolution, 
used  here  for  the  linear  solutions.  From  the  figure,  the  reason  for  the 
name  of  this  distribution  is  obvious  because  the  air  beam  forms  a small 
bump  on  the  tall  of  the  debris  distribution.  In  this  reference  frame, 
the  debris  Is  almost  stationary  so  the  air  velocity  Is  the  same  as  the 
relative  velocity.  This  velocity  Is  again  chosen  to  be  small  enough  so 
one-dimensional  analytical  theory  Is  applicable. 

The  analytical  solution  of  the  dispersion  relation  for  the  real 
part  of  the  frequency  by  the  power  series  method  mentioned  In  the 
previous  section  Is 
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In  Eqns  (103)  and  (104),  the  wave  phase  velocity  has  been  assumed  to  be 
much  greater  than  the  debris  thermal  speed.  Equation  (103)  Is  solved  by 
Iteration,  and  the  result  Is  Inserted  Into  Eqn  (104). 


The  test  calculations  are  done  using  the  plasma  distribution  shown 
In  Table  H-II.  Again,  the  air  density  Is  characteristic  of  that  found 
at  200  km.  The  analytical  solution  Is  found  for  k * 1 cm"1.  For  this 
wave  k,  only  stable  modes  exist.  The  analytical  and  computer  solutions 
for  one  such  mode  are 

to  » sec'1  (analytical)  (105) 


c*j  - S.iClo1^  — sec~'  (numerical)  (106) 

The  agreement  Is  good  with  the  discrepancy  being  due  to  the  warmth  of  the 
air  beam  In  the  computer  calculation.  The  computer  solution  gives  un- 
stable modes  for  wave  k's  from  0.075  cm"1  to  0.75  cm"1  with  the  fastest 
growing  mode  at  0.5  cm"1.  The  phase  velocities  of  the  growing  modes  are 
in  the  positive  slope  resonant  region  shown  in  Fig.  7.  Note  the  real 
part  of  the  frequency  Is  just  below  the  Ion  (debris)  plasma  frequency. 

Equal  Beam  with  Magnetized  Electrons.  The  final  linear  test  case 
is  similar  to  the  first  except  the  electrons  are  cold  and  magnetized. 

The  analytical  solution  uses  the  cold  hydrodynamic  result,  Eqn  (93), 
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modified  to  account  for  the  electron  magnetization.  Because  the  ions 
are  unmagnetized,  the  only  change  Is  to  the  electron  dielectric.  There- 
fore, when  Eqns  (53)  and  (76)  are  compared,  Equation  (93)  becomes 


ek  « [ — L 

r l (~  + 


The  electron  dielectric  also  results  directly  from  Eqn  (20)  when  Eqns 
(22),  (26),  (27),  and  (29)  are  used.  Equation  (107)  can  be  rewritten 
using  the  definition  of  the  lower  hybrid  frequency,  Eqn  (69): 


(107) 


I * ^ 


The  solution  of  Eqn  (108)  follows  that  of  Eqn  (93): 


>LH  [(o+Avf  + 


(108) 


= X v + - 1 


+ H CO 


(109) 


Analogous  to  Eqn  (99),  the  unstable  wave  k range  is 


O £ X < 


(110) 


The  test  plasma  is  shown  in  Table  H-III.  For  this  case,  the 


unstable  k range  is 


o < 


X < (.c 


and  the  complex  frequency  from  Eqn  (109)  for  k = 4.8  cm-1  is 
60  - 0.0  XiT//)  Sec*' 

The  computer  solutions  corresponding  to  Eqns  (111)  and  (112)  are 

3.r  < X < C~r 

u>  - 6.%(  icf)  sec*' 


(111) 


(112) 


(113) 


(114) 


The  agreement,  again.  Is  quite  good  with  the  discrepancy  in  the  lower 
end  of  the  wave  k range  being  due  to  the  resonant  effect  described 
earl  ler. 
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The  computer  solutions  for  this  test  case  have  two  very  Interesting 
features.  First,  the  unstable  modes  are  all  found  in  the  extremely  small 
angular  range:  5 degrees.  The  cold  electrons  freely  move  along 

the  magnetic  field  stabilizing  the  waves  with  even  small  components  of 
motion  parallel  to  the  field.  And  second,  the  wave  oscillation  frequen- 
cies quickly  Increase  as  the  propagation  direction  moves  away  from  being 
exactly  perpendicular  to  the  field.  As  $ increases,  the  wave  type  changes 
from  an  Ion-Ion  wave  to  a modified  two-stream  wave.  The  relationship 
between  the  Ion-Ion  and  modified  two-stream  wave  types  is  described  later 
in  this  chapter. 

Test  of  the  Evolution 

The  test  of  the  computer  evolution  of  the  ion  distributions  Is  done 
for  the  bump-on-the-tall  plasma  distribution  described  earlier.  Only  an 
unmagnetized  test  case  needs  to  be  done  because  the  Ion  diffusion 
equations  are  unchanged  when  the  electrons  become  magnetized.  In  fact, 
the  only  place  the  effect  of  electron  magnetization  appears  In  the 
quasi linear  equations  Is  In  the  electron  dielectric  term  of  the  disper- 
sion relation.  (The  magnetized  dispersion  relation  was  tested  In  the 
previous  section.)  In  the  test  evolution,  the  relative  velocity  Is 
chosen  to  be  far  enough  below  the  debris  acoustic  velocity  so  the  fastest 
growing  wave  modes  are  In  the  streaming  direction.  Thus,  the  computer 
evolution  should  be  approximately  one  dimensional. 

A one-dimensional  analytical  solution  of  the  quasll Inear  equations 
follows  Davidson  (Ref  14:  174-182)  and  is  In  Appendix  I.  This  solution 
gives  the  time  asymptotic  behavior  of  the  Ion  distributions  after  quasl- 
linear  saturation  has  occurred.  A plot  showing  qualitatively  the  nature 
of  the  asymptotic  solution  Is  Included  In  Fig.  7.  The  figure  shows  the 
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formation  of  a plateau  in  the  resonant  region  and  a slight  heating  of 
the  debris  ions.  No  further  changes  in  the  distributions  from  the 
configuration  shown  is  possible  according  to  quasil inear  theory. 

The  test  evolution  is  done  using  the  distribution  in  Table  H-IV. 


Unstable  wave  modes  exist  for 

1,1  £ £ 3.7 

(115) 

with  the  fastest  growing  mode: 

u>  ■=.  * U 6.3  (lo6)  it c~' 

(116) 

at  -A.  yg  = 3, 3.  < ***” 1 

(117) 

sfx  j r o.O  c »*i  ~ 1 

(118) 

A complete  listing  of  all  of  the  wave  modes  used  in  tl  test  evolution 
is  in  Table  J-I  In  Appendix  J.  Tabulated  in  this  appendix  are  all  of 
the  wave  modes  used  to  begin  each  of  the  dissertation  evolution  calcula- 
tions. These  modes  are  derived  from  solutions  of  the  linear  dispersion 
relation  using  the  initial  plasma  distributions.' 

The  possible  appearances  of  new  unstable  wave  modes  late  in  the 
evolution  calculations  are  Ignored  in  this  and  all  other  evolutions  In 
the  dissertation.  The  initial  modes  must  accumulate  spectral  energies 
several  orders  of  magnitude  above  the  initial  thermal  values  before  the 
distributions  begin  to  change.  This  large  head  start  Insures  the  In- 
significance of  any  late  appearing  modes  In  modifying  the  plasma.  Even 
among  the  initially  growing  modes,  the  fastest  growing  mode  quickly 
dominates  by  accumulating  a spectral  energy  much  larger  than  all  of  the 
rest. 

As  can  be  seen  from  Table  J-I,  all  of  the  growing  modes  propagate 
within  40  degrees  of  the  streaming  direction  (the  x direction)  so  the 
computer  quasilinear  solution  starting  with  these  modes  is  essentially 
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one  dimensional.  The  phase  velocities  over  cos^  of  these  modes  lie  In 
the  range: 

a.5"OoT)  1 ~ 2. (to1}  c (119) 

<-c«  Cr 

which  corresponds  to  the  positive  slope  of  the  combined  distribution 
shown  in  Fig.  7.  This  range  Is  In  the  resonant  region  which  changes  the 
most  rapidly  during  the  quasi! Inear  evolution. 

The  evolution  results  are  shown  In  Figs.  8 and  9.  In  these  figures 
and  In  all  such  figures;  the  debris  distribution,  the  air  distribution, 
and  their  sum  are  plotted  In  a three-dimensional  perspective  view.  The 
electrons  are  not  plotted  as  they  always  remain  In  their  initial  Max- 
wellian distribution.  The  plane  is  the  positive  uz  half  of  the  square 
which  represents  the  positive  vz  half  of  the  velocity  plane  as  is 
described  in  Appendix  F.  (The  plane  in  the  figure  is  the  computer 
velocity  mesh.)  The  point  of  view  is  above  the  positive  vz  axis  looking 
downward  and  in  the  negative  vz  direction.  The  front  and  back  edges  of 
the  plane  represent  vz  = and  vz  = 0 respectively.  The  left  and  right 
edges  represent  vx  = and  vx  = -<**  respectively.  The  heights  of  the 
distributions  are  normalized  to  the  original  Maxwellian  heights.  This 
normalization  is  held  constant  throughout  an  evolution  so  a change  in 
apparent  height  is  real.  Because  in  this  type  of  plot  appearance  of  the 
Initial  Individual  debris  and  air  Maxwell ians  Is  the  same  for  any  type 
of  plasma,  the  Initial  ion  distribution  plots  corresponding  to  Fig.  8a 
and  8b  are  omitted  for  the  remainder  of  the  dissertation.  The  computer 
program  constructs  these  three  types  of  plots  and  lists  the  distribution 
values  over  the  velocity  mesh  at  the  end  of  each  phase  of  the  evolution. 
These  plots  and  listings,  which  are  available  from  the  author,  are  used 


Fig.  8.  Initial  Distributions  for  the  Burr.p-on-the-Tail  Test 
Evolution. 
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to  check  on  the  general  progress  of  the  evolutions.  Because  the  plots 
are  Intended  to  show  the  qualitative  changes  to  the  Ion  distributions 
during  the  evolutions,  the  numerical  values  of  the  plotted  points  are  not 
Important,  and  hence,  are  omitted. 

In  Figs.  8 and  9,  the  Initial  distributions  and  the  final  distribu- 
tions calculated  are  respectively  shown.  The  final  distributions  are 
those  at  normalized  time  = 120  where  the  normalized  time  Is  In  units  of 
the  debris  plasma  oscillation  time, 

X,.,„  = "pV  (120) 

(This  time  normalization  Is  used  frequently  throughout  the  remainder  of 
the  dissertation.)  At  that  time,  the  wave  modes  were  rapidly  changing 
from  growing  modes  into  damped  modes,  and  the  anomalous  transport 
coefficients  were  declining.  Therefore,  Figure  9 shows  the  distributions 
close  to  saturation. 

The  computer  evolution  results  shown  In  Fig.  9 agree  well  both  with 
the  analytical  solution  results  shown  in  Fig.  7 and  with  the  assumptions 
used  In  the  dissertation  calculations.  A plateau  does  form  In  the  reso- 
nant region.  The  temperature  of  the  debris  increases  0.4  percent  which 
approximately  agrees  with  the  result,  0.2  percent,  from  Eqns  (1-32)  and 
(1-34).  This  change  In  the  nonresonant  particle  energy  Is  of  the  order 
of  the  wave  energy  as  predicted  by  Eqn  (1-36).  (The  discrepancy  Is 
mostly  due  to  the  large  amount  of  heating  In  the  tall  of  the  debris 
distribution.  The  tall  is  In  the  resonant  region,  but  the  computer 
temperature  calculation  Includes  the  tail  with  the  nonresonant  particles. 
In  fact,  the  change  In  the  resonant  particle  energy  Is  not  computed 
because  It  Involves  velocity  moments  of  parts  of  the  debris  and  air 
distributions.) 
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The  evolution  is  complete  In  about  10"®  seconds  after  the  ion  beams 
have  moved  40  cm  relative  to  each  other.  This  distance  (the  coupling 
length)  is  much  smaller  than  a typical  debris  bubble  mixing  region  size 
of  a few  kilometers  so  the  homogeneous  approximation  Is  a good  one.  The 
evolution  time  Is  much  less  than  the  10+*  seconds  needed  for  collisions 
to  alter  the  distributions  significantly  (Ref  28:  133),  yet  it  Is  greater 
than  the  time  from  Eqn  (A-31)  needed  for  the  neglected  free  streaming 
terms  to  damp  away. 

The  total  wave  energy  throughout  the  entire  evolution  Is  less  than 
the  debris  and  air  thermal  energies  so  Eqn  (45)  holds.  (Incidentally, 
most  of  the  total  particle  energy  is  in  the  electrons.)  Therefore,  the 
assumption  that  the  electric  fields  are  first  order  quantities  so  the 
higher  order  mode  coupling  terms  In  Eqn  (A-17)  can  be  ignored  is 
reasonable.  The  wave  energy  normalized  to  the  initial  total  ion  energy 
is  plotted  versus  time  in  Fig.  10.  As  can  be  seen  from  the  figure, 
most  of  the  evolution  consists  of  an  approximately  linear  phase  where 
the  wave  energy  grows  exponentially.  Late  in  the  evolution,  the  cal- 
culation becomes  nonlinear  as  the  changing  distributions  slow  the  wave 
energy  growth  rate.  This  behavior  is  typical  of  the  general  nature  of 
all  of  the  evolutions  where  nothing  much  happens  for  most  of  the  time, 
and  then  suddenly,  a lot  begins  to  happen  very  quickly. 

Because  late  in  the  evolution  the  total  wave  energy  approaches  the 
air  thermal  energy,  the  possibility  of  particle  trapping  should  be  in- 
vestigated. Particle  trapping  occurs  when  an  Ion  Is  trapped  by  a single 
wave  and  is  carried  along  with  it.  Trapping  is  a highly  nonlinear  process 
so  the  quasi 11  near  equations  break  down  when  it  occurs.  The  time  needed 
for  a singly  ionized  air  ion  to  be  trapped  is  (Ref  14;  54-56) 
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TOTAL  TIME  X DEBRIS  PLASMA  FREQUENCT 


Fig.  10.  Normalized  Total  Wave  Energy  Density  for  the  Bump-on- the 
Tall  Test  Evolution. 


T 

fraf/,,r 3 


(121) 


J ( lo~iS)  jfc 


(122) 


where  Ww  Is  the  average  electric  field  energy  density  In  the  few 
fastest  growing  wave  modes.  (Most  of  the  total  wave  energy  Is  In  these 
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modes.)  The  field  energy  Is  taken  at  Its  maximum  value  (the  worst  case). 
The  time  in  Eqn  (122)  Is  about  the  same  as  the  evolution  completion  time 
so  particle  trapping  cannot  begin  to  occur  until  the  evolution  Is  essen- 
tially complete. 

The  maximum  momentum  transport  coefficient  Is  0.17  This 

coefficient  corresponds  to  the  effective  anomalous  collision  frequency: 

V£  - **  (123) 

'“'A  < VA  ~ 

X 'l  (/oO  sec'1  (124) 

where  the  numerator  Is  obtained  from  Eqn  (36)  and  the  computer  results. 
This  collision  frequency  Is  about  0.2  percent  of  the  debris  plasma  fre- 
quency. (Note  throughout  the  remainder  of  this  chapter,  the  anomalous 
transport  coefficient  results  will  always  be  quoted  In  terms  of  an 
effective  momentum  transfer  collision  frequency  for  air  normalized  to 
the  debris  plasma  frequency.)  The  collision  time,  which  Is  the  recipro- 
cal of  Eqn  (124),  Is  much  less  than  the  colllslonal  relaxation  time 
mentioned  above.  Therefore,  the  instability  modifies  the  plasma  much 
faster  than  binary  collisions  in  agreement  with  the  collisionless  assump- 
tion of  Vlasov  theory. 

The  temporal  behavior  of  the  normalized  anomalous  collision  fre- 
quency shown  In  Fig.  11.  The  exponential  growth  early  In  the  evolution 
Is  clear  as  well  as  the  fall  off  late  in  the  evolution  as  the  distri- 
butions begin  to  saturate.  The  fall  off  is  due  to  the  rapid  decline  of 
both  the  mode  growth  rates  and  Imaginary  parts  of  their  Ion  dielectrics. 
The  average  air  velocity  normalized  to  the  debris  acoustic  speed  Is  also 
shown  In  the  figure.  The  decrease  In  the  air  velocity  Is  the  main  energy 
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TOTAL  TIME  X OEBRIS  PLASMA  FREQUENCY 

Fig.  11.  Anomalous  Collision  Frequency  and  Average  Air  Velocity 
for  the  Bump-on- the-Tail  Test  Evolution. 

source  driving  the  Instability. 

The  temporal  behavior  of  the  initially  fastest  growing  mode  is 
shown  In  Fig.  12.  The  real  and  Imaginary  parts  of  Its  frequency  normal- 
ized to  the  debris  plasma  frequency  are  shown.  Again,  the  approximately 
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Fig.  12.  Temporal  Behavior  of  the  Initially  Fastest  Growing  Mode 
for  the  Bump-on- the-Tail  Test  Evolution. 

linear  phase  followed  by  the  nonlinear  phase  is  evident.  The  growth  rate 

falls  because  the  changes  in  the  distributions  tend  toward  stabilization. 

The  decrease  in  the  real  part  reflects  the  decrease  In  the  resonant 

region  velocity  as  the  air  heats  and  slows  down.  The  real  part  of  the 

frequency  decreases  to  maintain  the  resonance  given  by  Eqn  (18).  The 

small  oscillation  in  this  curve  is  caused  by  a small  inflection  In  the 

air  distribution  slope  passing  the  velocity  resonant  with  this  mode. 

The  sudden  change  In  the  slope  of  the  curve  for  the  imaginary  part 

near  time  = 110  in  Fig.  12  Is  due  to  the  calculation  beginning  to  break 

down.  The  wave  modes  very  late  in  the  evolution  were  not  updated  often 
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enough.  The  development  of  a method  of  determining  the  frequency  of 
mode  update  during  an  evolution  Is  not  a trivial  task.  Several  Itera- 
tions were  required  to  get  a suitable  method.  The  method  must  update 
infrequently  early  In  the  evolution  when  the  energies  in  the  wave  modes 
are  small  and  update  frequently  late  In  the  evolution  when  the  distribu- 
tions are  changing  rapidly.  The  distributions  In  a calculation  which 
breaks  down  develop  oscillations  essentially  between  plus  and  minus 
Infinity  within  a few  hundreths  of  an  e- folding  time.  The  update  control 
used  in  the  trial  evolution  was  an  earlier  attempt  at  the  control  algo- 
rithm. Fortunately,  a small  degree  of  breakdown  has  little  effect  in 
the  maximum  anomalous  transport  coefficients. 

In  conclusion,  the  computer  evolution  appears  physically  sound  and 
agrees  with  an  analytical  solution.  The  results  also  seem  to  be  rela- 
tively accurate  and  consistent.  The  anomalous  momentum  transfer 
coefficients  for  the  debris  and  air  should  be  equal  and  opposite.  They 
are  to  two  significant  figures.  The  electron  transport  coefficients  are 
much  smaller  than  those  for  the  Ions  so  the  adiabatic  electron  assump- 
tion is  a good  one.  The  numerical  version  of  the  dispersion  relation 
applied  to  the  initial  Maxwellian  distributions  gives  results  which 
agree  with  the  analytical  dispersion  relation  results.  The  oscillation 
frequencies  are  the  same  to  two  or  three  places,  and  the  growth  rates 
usually  agree  to  about  one  or  two  places.  Occasionally,  however, 
discrepancies  of  up  to  a factor  of  two  appear  in  the  growth  rates. 
Particle,  momentum,  and  energy  densities  checked  by  CONSERV  are  conser- 
ved to  within  one  percent.  The  fast  integration  routine  for  the  diffu- 
sion coefficients,  COEFINT,  Is  not  very  accurate  compared  to  SINGRT. 

Yet,  the  net  diffusion  Is  the  same  with  the  maximum  anomalous  transport 
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coefficients  from  an  evolution  using  each  routine  agreeing  to  two  places. 
Evidently,  because  so  many  diffusion  coefficient  Integrals  are  done,  the 
Inaccurate  COEFINT  results  average  out  to  those  of  SINGRT.  In  general, 
the  computer  evolution  results  seem  to  be  accurate  to  at  least  within  a 
factor  of  two. 

Debris-Air  Evolution  Resul ts 

The  debris-air  evolutions  are  divided  into  two  categories  according 
to  whether  or  not  the  electrons  are  magnetized.  The  reasons  for  these 
two  categories  and  the  criteria  for  placing  a given  debris-air  problem 
Into  one  of  them  are  covered  In  Chapters  I and  II.  The  evolutions  using 
air  densities  typical  of  altitudes  of  200  km  fall  Into  the  unmagnetized 
category  because  the  resulting  electron  distribution  can  be  expected  to 
overpo wer  any  magnetic  fields  which  are  present.  For  air  densities 
typical  of  600  km,  the  reverse  holds  so  such  evolutions  fall  into  the 
magnetized  category.  Also,  at  600  km  the  electrostatic  criterion,  Eqn 
(78),  holds  quite  well  for  reasonable  magnetic  field  strengths.  The 
evolution  results  from  sample  problems  In  each  of  these  categories  are 
described  below. 

Unmagnetized  Evolution  Results.  The  first  sample  unmagnetized 
debris-air  problem  is  the  same  bump-on-the-tall  distribution  used  above 
except  with  a relative  velocity  greater  than  the  ion  acoustic  speed. 

The  species  distributions  and  relevant  parameters  are  shown  In  Table 
H-V.  This  distribution  Is  chosen  so  that  two-dimensional  effects  are 
emphasized  and  so  the  results  can  be  compared  to  and  contrasted  with  the 
largely  one-dimensional  test  case.  The  neglect  of  the  magnetic  field  Is 
a good  assumption  because  a reasonable  25-gauss  field  and  a typical 

A O 

electron  temperature  of  3(1(T)  eV  give  an  electron  of  2(10).  The 
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linear  solutions  for  the  unstable  modes  are  shown  in  Table  J-II. 

The  results  of  this  evolution  are  shown  in  Figs.  13  and  14.  The 
distributions  shown  are  respectively  the  Initial  one  for  the  ions  and 
those  at  the  time  of  the  maximum  momentum  transfer  rate.  The  fastest 
growing  mode  initially  propagates  at  a 73-degree  angle.  The  heating  of 
the  air  at  approximately  that  angle  to  the  negative  vx  direction  Is 
apparent  In  Fig.  14.  The  plateau  shown  In  Fig.  7 does  not  form.  When 
the  evolution  is  finished,  the  angle  of  maximum  growth  decreases  to  69 
degrees  but  the  large  accumulated  wave  energy  In  the  Initially  fastest 
growing  mode  still  dominates  the  evolution. 

The  maximum  value  of  the  momentum  transfer  rate  Is  0.078 
which  is  less  than  half  of  that  for  the  one-dimensional  evolution  even 
though  the  Initial  wave  growth  rates  are  nearly  the  same.  The  momentum 
transfer  rate  In  terms  of  the  normalized  collision  frequency  Is  plotted 
In  Fig.  15.  The  normalized  average  air  velocity  is  also  shown.  The 
dotted  line  extension  to  the  collision  frequency  curve  is  the  result  from 
an  earlier  attempt  at  the  evolution  which  began  to  break  down  just  before 
the  peak  was  reached.  This  result  clearly  shows  the  relative  insensi- 
tivity of  the  collision  frequency  on  the  exact  details  of  an  evolution. 
This  breakdown  Is  also  shown  In  the  next  figure.  Fig.  16.  This  figure 
also  reveals  that  the  real  part  of  the  wave  mode  frequency  increases 
with  time.  This  Is  true  for  most  of  the  evolutions,  but  the  reason  for 
this  Is  not  clear.  Of  course,  the  imaginary  part  always  decreases  with 
time. 

Because  the  unstable  waves  for  this  distribution  (where  VR/vDgc  = 
2.6)  propagate  nearly  perpendicular  to  the  streaming  direction,  the  dls- 
trlfc  tlon  temperature  In  the  vz  direction  Increases  dramatically.  The 
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Fig.  15.  Anomalous  Collision  Frequency  and  Average  Air  Velocity 
for  the  Large  Relative  Velocity  Bump-on- the-Tail 
Evolution. 


debris  and  air  v2  temperatures  by  the  time  the  maximum  momentum  transfer 


rate  is  reached  respectively  Increase  0.35  percent  and  110  percent.  The 


vx  temperatures  actually  decrease  0.02  percent  and  11  percent  respec- 


tively These  results  contrast  with  the  strong  vx  heating  in  the  one- 


dir.iensional  (Vj,/vCac  = 0.78)  evolution.  The  corresponding  values  to 


those  above  for  the  vz  temperature  increases  are  0.01  percent  and  -5.1 


percent,  and  those  for  v^  are  0.35  percent  and  91  percent.  The  one- 


dimensional anomalous  heating  rates  are  somewhat  less  than  those  for  the 
two-dimensional  evolution.  These  observations  indicate  that  at  the 


68 


I 


COLLISION  FREQUENCY/ DEBRIS  PLASMA  FREQUENCY 


AFIT/DS/PH/78-1 


>- 


TOTAL  TIME  X DEBRIS  PLASMA  FREQUENCY 

Fig.  16.  Temporal  Behavior  of  the  Initially  Fastest  Growing 
Mode  for  the  Large  Relative  Velocity  Bump-on- the- 
Tail  Evolution. 

higher  relative  velocities  more  of  the  wave  energy  heats  the  distribu- 
tions at  large  angles  to  the  streaming  direction  but  less  momentum 
transfer  occurs.  The  Importance  of  the  use  of  two  dimensionality  in  the 
evolutions  is  apparent. 

The  coupling  length  for  the  two-dimensional  evolution  is  130  cm 
which  is  still  much  smaller  than  the  debris-air  mixing  region.  The  ratio 
of  the  total  wave  energy  at  the  time  of  the  peak  momentum  transfer  rate 
to  the  initial  thermal  wave  energy  is  6.3(10  ).  The  evolution  took  50 
min  of  computer  time  while  the  one-dimensional  evolution  took  2 1/4  hrs. 
The  difference  is  mostly  due  to  the  relatively  large  number  of  wave 
modes  followed  in  the  one-dimensional  case. 
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The  next  two  evolutions  are  the  2C0-km  simulations  which  are  done 
for  low  and  high  relative  velocities.  The  low  velocity  distribution 
parameters  are  shown  in  Table  H-VI,  and  those  for  the  high  velocity  run 
are  In  Table  H-VI I . The  corresponding  linear  solutions  are  in  Tables 
J-III  and  J-IV.  Before  these  distributions  were  evolved,  the  possible 
existence  of  the  Ion  acoustic  Instability  between  the  electrons  and 
either  of  the  Ion  beams  was  examined  (Ref  19:  476-478).  The  ion  acoustic 
modes  do  not  exist,  however,  because  the  warm  temperatures  of  the  Ion 
beams  shut  off  the  instability.  (The  Ion  acoustic  instability  will 
occur  only  if  the  ions  are  cold  and  have  a large  velocity  relative  to 
the  electrons  so  the  ion  Landau  damping  does  not  overwhelm  the  weak 
Landau  growth  from  the  electrons.) 

The  results  of  these  evolutions  are  shown  in  Figs.  17  through  20. 
Again,  the  figures  show  the  initial  combined  ion  distributions  only  and 
the  debris,  air,  and  combined  ion  distributions  at  the  time  of  the 
maximum  momentum  transfer  rate.  The  results  in  the  figures  again  show 
the  one  dimensionality  of  the  low  relative  velocity  (VR/vDaC  * 1.2) 
evolution  and  the  two  dimensionality  of  the  high  velocity  (Vp/vpac  " 3.8) 
evolution.  The  low  velocity  results  in  Fig.  18  show  predominately  one- 
dimensional heating.  This  occurs  because  the  initially  fastest  growing 
mode  propagates  at  zero  degrees.  Both  the  debris  and  air  Ions  experience 
resonant  heating  with  the  flattened  areas  of  their  distributions  cor- 
responding to  the  phase  velocities  of  the  growing  modes.  The  one- 
dimensional heating  Is  reflected  in  the  changes  In  the  debris  and  air 
temperatures.  The  changes  In  the  vx  direction  are  8.6  percent  and  18 
percent,  and  the  changes  In  the  v2  direction  are  -0.45  percent  and  -2.7 
percent.  The  relative  velocity  decreases  1.5  percent  providing  the 
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energy  to  drive  the  instability.  The  high  velocity  results  in  Fig.  20 
show  heating  at  large  angles  to  the  streaming  direction  because  the 
fastest  growing  modes  propagate  approximately  between  55  and  70  degrees. 
The  flattened  area  on  the  air  distribution  corresponds  to  one  of  the 
resonant  regions  shown  in  Fig.  4.  No  resonant  heating  of  the  debris 
occurs.  The  changes  in  the  debris  and  air  temperatures  in  the  vx 
direction  are  1.3  percent  and  3.5  percent,  and  the  changes  in  the  vz 
direction  are  5.3  percent  and  12  percent.  The  relative  velocity 
decreases  0.39  percent. 

The  normalized  anomalous  collision  frequencies  for  both  of  the  200- 
km  simulations  are  graphed  in  Fig.  21.  The  maximum  anomalous  collision 
frequency  for  the  low  velocity  simulation  is  roughly  an  order  of  mag- 
nitude higher  than  that  for  the  high  velocity  run.  The  maximum  for  the 
low  velocity  also  occurs  in  about  one  half  of  the  time  of  the  other. 

This  result  is  also  reflected  in  the  larger  difference  in  coupling 
lengths  than  can  be  accounted  for  by  the  difference  in  relative  veloci- 
ties. The  lengths  for  the  low  and  high  velocity  runs  are  31  and  160  cm. 
The  anomalous  energy  transport  rates  for  the  high  velocity  run  are  about 
a factor  of  two  higher  than  the  low  velocity  results. 

The  imaginary  parts  of  the  Initially  fastest  growing  modes  in  the 
two  simulations  are  shown  In  Fig.  22.  Even  though  the  low  velocity  run 
ends  in  half  the  time  and  has  an  order  of  magnitude  larger  anomalous 
collision  frequency.  Its  initial  growth  rate  is  only  a third  larger. 

The  ratios  of  the  total  wave  energies  at  the  end  of  the  two  simulations 
to  those  at  their  beginnings  are  1.6(10^)  and  3.2(10^)  respectively  for 
the  low  and  high  velocity  runs.  The  respective  computer  running  times 

i 

were  2.4  and  1.3  hrs. 
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Fig.  21.  Anomalous  Collision  Frequencies  for  the  200-km 
Simulations. 

Magnetized  Evolution  Results.  The  magnetized  cases  simulating  the 
600-km  problem  are  somewhat  more  complicated  to  evolve  for  several 
reasons.  First,  the  modified  two-stream  Instability  often  competes  with 
the  Ion-Ion  Instability  (Refs  3:  53-56  and  29).  Fortunately,  only  the 
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Fig.  22.  Temporal  Behavior  of  the  Initially  Fastest  Growing 
Modes  for  the  20C-km  Simulations. 

ion-ion  instability  operates  for  the  high  electron  temperatures  used  in 
the  simulations.  Next,  the  electrostatic  criterion  also  has  to  be 
satisfied.  This  requires  some  juggling  of  the  relative  velocities,  the 
magnetic  field  strength,  and  the  Ion  densities  (the  simulation  altitude) 
to  insure  that  the  criterion  Is  rigorously  satisfied  for  the  evolutions 
used  in  the  dissertation.  Finally,  the  magnetized  evolutions  use  more 
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computer  time  because  the  dispersion  relation  is  more  complicated  and  the 
evolutions  take  longer  to  reach  the  peak  momentum  transfer  rate. 

The  mixture  of  modified  two-stream  and  ion-ion  modes  which  appears 
when  the  electron  temperature  is  reduced  is  shown  in  Table  J-V.  This 
table  shows  the  linear  solutions  for  the  distribution  shown  in  Table 
H-VIII  except  the  electron  temperature  in  the  vz  direction  is  an  order 
of  magnitude  lower  than  the  tabulated  value.  The  type  of  mode  is  deter- 
mined by  the  relative  magnitudes  of  the  imaginary  parts  of  the  dielec- 
trics. For  example,  the  mode  at  kx  = 0.1064  cm"1  and  kz  = 0.0  cm"1  is 
ion-ion  because  the  imaginary  part  of  the  electron  dielectric  is  much 
less  than  those  of  the  ion  dielectrics.  Conversely,  the  mode  at  kx  = 
0.1064  cm-1  and  kz  = 3.125(10“*)  cm"1  is  a modified  two-stream  mode 
between  the  air  and  the  electrons.  Note  also,  the  extremely  small  angu- 
lar range  of  the  unstable  modes. 

The  distributions  evolved  are  in  Tables  H-VIII  and  H-IX  respectively 
for  low  and  high  relative  velocities.  The  corresponding  linear  solutions 
are  in  Tables  J-VI  and  J-VII.  All  of  the  wave  modes  shown  are  clearly 
ion-ion.  Before  the  distributions  were  evolved,  the  possible  existence 
of  modified  two-stream  modes  propagating  to  the  "left”  (See  Fig.  6)  was 
also  checked.  No  such  modes  were  found  so  all  of  the  unstable  modes  for 
the  distribution  are  ion-ion,  and  therefore,  the  electrons  are  adiabatic. 

The  evolution  results  are  shown  in  Figs.  23  through  26.  The  low 
relative  velocity  (Vp/Vpa^  e 0.042)  evolution  results  are  very  similar 
to  those  of  the  low  velocity  evolution  for  200  km.  The  heating  is  one 
dimensional,  and  the  phase  velocities  of  the  unstable  modes  and  the 
resonant  areas  of  the  distributions  are  the  same  as  those  in  the  200-km 
low  velocity  case.  The  changes  in  the  debris  and  air  temperatures  in 
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Fig.  23.  Initial  Ion  Distribution  for  the  600-km  Simulation 
with  a Low  Relative  Velocity. 


/CTN  24a.  Debris 


bffgs 


Fig.  24.  PeakVc  Distributions  for  the  600-km  Simulation  with 
a Low  Relative  Velocity. 
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the  direction  are  20  percent  and  35  percent,  and  the  changes  In  the 
v2  direction  are  -1.8  percent  and  -5.5  percent.  These  changes  are  much 
larger  than  the  200-km  results.  The  relative  velocity  decreases  3.3 
percent. 

The  high  velocity  results  are  quite  different  from  the  corresponding 
200-km  case.  The  relative  velocity  used  is  between  the  debris  Alfven 
and  acoustic  speeds  (vp/vDal^  =0.15  and  VR/vDac  = 4.5)  so  one-dimen- 
sional magnetized  modes  and  two-dimensional  unmagnetized  modes  both  exist. 
This  division  Is  shown  in  Table  0-VII,  the  linear  solutions.  The  mag- 
netized modes  propagate  In  the  streaming  direction,  and  the  unmagnetized 
modes  propagate  In  an  angular  range  near  60  degrees.  The  unmagnetlzed 
inodes  are  very  similar  to  the  wave  modes  for  the  high  velocity,  200-km 
case.  Because  the  magnetized  modes  have  the  larger  growth  rates,  the 
heating  is  predominately  In  the  Vx  direction  as  shown  In  Fig.  26.  The 
changes  In  the  debris  and  air  temperatures  In  the  v direction  are  110 
percent  and  190  percent,  and  the  changes  In  the  vz  direction  are  -18 
percent  and  -36  percent.  The  heating  of  the  distributions  Is  fairly 
uniform  with  little  resonant  heating  being  observed.  The  relative 
velocity  decreases  2.7  percent. 

The  normalized  anomalous  collision  frequencies  for  these  simulations 
are  shown  In  Fig.  27.  The  peak  collision  frequency  for  the  high  velocity 
run  Is  about  four  times  as  large  as  that  for  the  low  velocity  run.  The 
next  figure.  Fig.  28,  shows  the  growth  rates  of  the  Initially  fastest 
growing  modes.  The  Initial  growth  rate  for  the  high  velocity  distribu- 
tion Is  twice  that  of  the  other.  The  ratios  of  the  total  wave  energies 

5 

at  the  end  of  the  simulations  to  those  at  their  beginnings  are  5.7(10  ) 
and  4.0(107). 
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Fig.  27.  Anomalous  Collision  Frequencies  for  the  600-km 
Simulations. 

The  coupling  lengths  for  the  low  and  high  velocity  runs  are  880  and 
1900  cm.  These  lengths  are  much  larger  than  an  Ion  gyro  radius  so  the 
assumption  that  the  Ions  are  field  free  breaks  down.  In  the  debris-air 
coupling  problem,  this  breakdown  means  that  a mechanism  based  on  curving 
Ion  paths  causing  the  coupling  (Ref  30)  is  more  likely. 

The  computer  running  times  for  these  evolutions  were  2.7  and  3.8 
hrs.  The  long  times  were  partly  due  to  the  use  of  an  earlier  Inefficient 
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Fig.  28.  Temporal  Behavior  of  the  Initially  Fastest  Growing 
Modes  for  the  600- km  Simulations. 

control  on  the  rode  updates. 

The  above  discussion  on  the  600-km  simulation  results  completes  the 
presentation  of  the  dissertation  evolution  calculations.  In  the  first 
few  sections  of  this  chapter,  the  testing  of  the  computer  program  written 
to  do  the  evolutions  was  described.  Included  in  that  description  was  an 
evaluation  as  to  how  well  the  computer  results  satisfy  the  assumptions 
of  quasil Inear  theory.  In  the  remaining  sections  of  the  chapter,  the 
results  of  the  simulation  calculations  were  presented.  These  results 
are  summarized  and  evaluated  in  the  next  chapter. 
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V.  Discussion  and  Recommendations 

First  in  this  chapter,  the  results  of  the  evolutions  done  for  the 
dissertation  are  surunarlzed.  Next,  the  results  are  checked  for  consis- 
tency both  with  the  assumptions  made  In  Chapter  I and  with  the  assump- 
tions of  quaslllnear  theory.  And  finally  In  the  chapter,  some  recom- 
mendations for  further  study  are  made. 

Summary  of  Evolution  Results 

Several  general  features  of  the  evolutions  are  apparent.  First, 
the  nature  of  the  evolutions  changes  dramatically  as  time  progresses. 

At  early  times,  nothing  much  happens  because  the  wave  energies  are  still 
small.  At  late  times,  the  Ion  distributions  begin  to  change  rapidly,  and 
the  growing  wave  modes  rapidly  begin  to  drop  out  becoming  damped  modes. 
The  bulk  of  the  computer  time  is  spent  In  the  latter  phases  of  the  evolu- 
tions. 

The  next  general  feature  of  the  evolutions  is  that  the  maximum 
anomalous  transport  rates  occur  before  saturation.  Although  these  rates 
are  proportional  to  the  wave  energy,  they  also  depend  on  the  distribution 
slopes  through  the  wave  growth  rates  and  plasma  dielectrics.  (See  Eqns 
(36)  and  (37).)  These  slopes  decrease  rapidly  late  In  the  evolution 
calculations  causing  the  transport  rates  to  peak  and  then  decline. 

The  evolutions  verify  the  characteristic  frequencies  and  velocities 
derived  In  Chapter  II.  The  unmagnetized  waves  do  tend  to  oscillate  near 
the  Ion  plasma  frequency.  This  tendency  Is  apparent  when  one  compares 
the  Initial  wave  oscillation  frequencies  for  the  bump-on- the- tall  cases 
In  Tables  0-1  and  0- II  with  the  debris  plasma  frequency  in  Table  H-IV. 
This  tendency  Is  less  apparent  when  the  Ion  beams  are  of  comparable 
densities  because  the  Interactions  of  the  debris-electron  oscillations 
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with  the  air-electron  oscillations  reduce  the  resulting  wave  frequencies. 
In  fact,  when  the  Ion  beams  are  equal  and  opposite  (o^p  * <o^),  tfle 
wave  frequencies  go  to  zero.  Similarly,  the  magetized  waves  do  tend  to 
oscillate  at  the  lower  hybrid  frequency.  Note  these  results  confirm 
the  reasonableness  of  the  time  normalization  In  Eqn  (120).  The  charac- 
teristic velocity  for  the  unmagnetized  waves  Is  the  ion  acoustic  speed 
because  the  wave  phase  velocities  do  tend  to  be  near  this  speed.  The 
actual  phase  velocities  are  somewhat  less  than  the  ion  acoustic  speed, 
however,  due  to  the  wave  Interaction  effect  described  above.  (Also,  see 
Eqn  (18).)  Eecause  the  waves  cannot  move  faster  than  the  Ion  acoustic 
speed,  the  unstable  wave  modes  for  plasma  distributions  with  high 
relative  velocities  must  propagate  at  large  angles  to  the  streaming 
direction.  Their  phase  planes  can  then  keep  up  with  the  resonant  areas 
as  shown  In  Fig.  4.  The  characteristic  velocity  for  the  magnetized  waves 
seems  to  be  the  Alfven  speed.  This  is  most  apparent  in  the  high 
velocity,  600-km  evolution.  In  that  evolution,  the  Initial  magnetized 
waves  propagate  faster  than  the  ion  acoustic  speed,  but  they  still  are 
In  the  streaming  direction  because  the  relative  velocity  Is  less  than 
the  Alfven  speed.  In  general,  the  relationship  between  the  character- 
istic velocity  and  relative  velocity  of  a given  plasma  type  determines 
. whether  the  evolutions  will  be  one  dimensional  or  two  dimensional. 

The  next  feature  of  the  evolutions  is  that  they  are  definitely  Ion- 
ion.  The  transfer  of  momentum  and  energy  from  the  Ions  to  the  electrons 
Is  negligible  compared  to  the  transfers  between  the  two  Ion  species. 

Another  feature  Is  that  the  peak  transport  rates  are  attained  In 
about  ten  e-foldlng  times  where  an  e-folding  time  Is  the  reciprocal  of 
the  largest  Initial  wave  growth  rate.  This  result  Is  quite  different 
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from  the  usual  "back-of-the-envelope"  time  to  quaslllnear  saturation  of 
one  e-foldlng  time.  This  time  to  attain  the  peak  transport  rates  also 
corresponds  to  the  time  when  the  wave  energy  begins  to  become  comparable 
to  the  Ion  thermal  energy. 

The  evolutions  are  characterized  by  strong  diffusion  In  the 
regions  of  the  distributions  resonant  with  the  waves  and  by  weak  dif- 
fusion In  the  nonresonant  regions.  The  flattened  areas  on  the  distribu- 
tions In  the  figures  in  the  previous  chapter  are  the  result  of  the  strong 
diffusion.  The  weak  diffusion  cannot  be  seen  In  the  figures,  but  It  Is 
apparent  In  the  numerical  printouts  of  the  distributions.  In  general, 
the  effect  of  the  nonresonant  diffusion  is  a slight  heating  of  the  Ion 
distributions  similar  to  that  shown  In  Fig.  7. 

The  evolutions  appear  to  be  numerically  stable.  Occasionally,  a 
small  dip  at  a single  mesh  point  in  an  ion  distribution  suddenly 
appeared  during  an  evolution  calculation  probably  due  to  an  error  In  the 
GEARB  solution.  This  dip  quickly  filled  In  during  the  next  few  phases 
of  the  evolution  thereby  proving  the  Inherent  stability  of  the  solution. 
Yet,  If  the  wave  modes  were  not  updated  often  enough,  an  evolution  would 
quickly  become  unstable  with  rapidly  growing  oscillations  in  the  particle 
distributions.  Clearly,  the  most  important  factor  on  numerical  stability 
is  the  frequency  of  mode  update. 

The  final  general  feature  of  the  dissertation  results  Is  that  more 
momentum  transfer  and  less  particle  heating  occurs  In  the  low  relative 
velocity  evolutions  than  in  the  high  velocity  ones.  This  phenomenon  Is 
reasonable  because  the  predominately  z direction  diffusion  In  the  high 
velocity  runs  is  symmetric  as  shown  in  Fig.  4.  A momentum  transfer  In 
the  positive  z direction  In  the  upper  resonant  region  Is  canceled  out 
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by  a momentum  transfer  In  the  negative  z direction  In  the  lower  region. 

The  net  effect,  therefore.  Is  a spreading  of  the  distribution  which 
Increases  the  z temperature.  The  low  velocity  diffusion,  on  the  other 
hand.  Is  predominately  In  the  x direction  as  shown  In  Fig.  3.  This 
diffusion  Is  not  symmetric  so  a net  transfer  of  momentum  does  occur.  The 
phenomenon  can  be  seen  In  Table  I which  Is  a summary  of  the  dissertation 
evolution  results.  The  tabular  entries  are  in  a suitable  form  for 
Insertion  Into  a CGS  set  of  hydrodynamic  equations.  The  Implication  of 
these  results  for  the  exospheric  nuclear  explosion  coupling  problem  Is 
that  the  icn-lon  Instability  is  not  effective  in  coupling  the  air  to  the 
debris  when  their  relative  velocity  is  high. 

Validity  of  Approximations 

The  assumptions  used  in  the  dissertation  quasillnear  calculations 
hold  fairly  well  In  most  cases.  The  electric  field  energy  in  the  unstable 
waves  does  remain  smaller  than  the  particle  thermal  energy  in  accordance 
with  the  weak  turbulence  criterion,  Eqn  (45).  The  wave  electric  fields, 
therefore,  are  of  first  order,  and  the  perturbation  approach  used  to 
derive  the  quasilinear  equations  is  valid.  The  electrostatic  approxima- 
tion criterion,  Eqn  (78),  holds  for  the  magnetized  distributions  used. 

(The  criterion  Is  not  pertinent  for  the  unmagnetized  cases  because 
coupling  between  the  electrostatic  waves  and  possible  electromagnetic 
waves  does  not  occur.)  The  wave  oscillation  frequencies  for  the 
magnetized  cases  do  fall  In  the  range  defined  by  Eqn  (29).  The  average 
electron  gyro  radii  are  small  so  Equation  (31)  is  satisfied.  The  ion 
radii,  however,  are  large  compared  to  the  mode  wavelength  as  required 
In  part  of  Eqn  (30),  but  they  are  small  compared  to  the  coupling  length 
violating  the  other  part  of  that  equation.  This  violation  means  the 
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magnetized  ion-ion  instability  does  not  reach  saturation  before  the 
plasma  is  altered  by  the  Ion  gyro  motions  around  the  magnetic  field 
lines.  The  coupling  lengths  and  unstable  mode  wavelengths  are  always 
less  than  a few  meters  so  the  homogeneous  approximation  Is  a good  one. 
(Recall  the  mixing  region  thickness  Is  of  the  order  of  a few  kilometers.) 
The  anomalous  collision  frequencies  are  always  much  larger  than  the 
binary  collision  frequencies  so  the  Vlasov  approach  is  valid.  The  phase 
mixing  times  for  the  free-streaming  terms  are  always  much  less  than  the 
evolution  times  (the  times  to  reach  the  maximum  transport  rates)  so  the 
free-streaming  terms  can  be  neglected.  The  trapping  times  are  of  the 
order  of  the  evolution  times,  however. 

In  general,  the  unmagnetized  quasilinear  evolutions  done  for  the 
dissertation  seem  to  be  valid  up  to  about  the  time  the  transport  coef- 
ficients peak  and  then  begin  to  decline.  After  that  time,  the  wave 
energy  begins  to  become  too  large,  particle  trapping  begins  to  occur, 
and  the  mode  growth  rates  begin  to  change  rapidly.  These  factors 
combine  to  make  evolution  calculations  at  times  much  beyond  the  peaks 

invalid  and  prohibitively  expensive.  These  problems  result  from  the 

fundamental  limitation  of  quasilinear  theory:  It  is  a weak  turbulence 

theory,  and  therefore,  it  cannot  handle  wave-wave  or  strong  wave-particle 

interactions.  Yet,  by  the  time  the  evolution  calculations  become  in- 
valid, their  main  objective,  which  is  a reasonably  accurate  set  of 
maximum  anomalous  transport  coefficients,  has  been  attained. 

The  magnetized  evolutions  appear  to  be  limited  by  the  ion  gyration 
times  (Jd^1)  which  are  much  less  than  the  evolution  times.  A simple 
solution  to  this  problem  is  not  readily  apparent  within  the  scope  of 
the  theory  in  the  dissertation.  Yet,  in  these  cases  too,  quasilinear 
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theory  has  served  a useful  purpose  because  the  calculated  transport 
coefficients  are  reasonable  estimates  of  the  actual  ones. 

Recommendations  for  Further  Study 

Two  of  the  more  obvious  immediate  extensions  of  the  dissertation 
effort  are  (1)  to  run  more  evolutions  to  determine  how  the  results  scale 
and  (2)  to  generalize  the  computer  program  so  it  can  do  instability 
evolutions  with  nonadiabatic  electrons.  A few  aspects  of  scaling  are 
fairly  straightforward  and  would  require  only  an  evolution  or  two  to 
verify.  For  example,  if  the  temperatures  of  the  species,  the  Ion 
relative  velocity,  and  the  ratio  of  the  debris  plasma  frequency  to  that 
of  the  air  are  all  held  constant  while  the  debris  density  varies;  the 
change  In  wave  mode  growth  rates  will  match  the  change  in  debris  plasma 
frequency.  Also,  the  spectrum  of  unstable  modes  will  always  propagate  in 
the  same  angular  range.  This  behavior  can  already  be  seen  in  the  un- 
magnetized modes  for  the  high  relative  velocity  200-km  and  600-km  runs 
in  Tables  J-IV  and  J-VII.  The  ratios  of  the  growth  rates  of  correspon- 
ding modes  are  all  the  same  and  are  equal  to  the  ratio  of  the  debris 
plasma  frequencies.  Therefore,  the  maximum  transport  cofficients  should 
scale  as  (See  Eqn  (11).) 
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The  computer  program  has  already  been  generalized  enough  so  It  can 
handle  unmagnetized  electron-ion  instabilities.  In  this  generalization, 
the  electron  dielectric  is  computed  numerically  using  the  same  procedures 
described  in  Appendices  F and  G for  the  always  unmagnetized  ions.  This 
numerical  computation  is  necessary  because  the  electrons  are  not  always 
adiabatic  in  electron-ion  evolutions.  More  importantly,  however,  the 
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program  needs  to  be  generalized  enough  so  it  can  evolve  the  modified 
two-stream  instability  because  the  wave  k's  for  this  instability  lie 
close  to  those  for  the  ion-ion  instability.  (See  Table  J-V.)  To  do 

] 

this,  the  routines  which  calculate  the  magnetized  electron  dielectric 
have  to  be  modified  to  accommodate  a non  Maxwellian  electron  distribution. 

Another  relatively  straightforward  extension  of  the  dissertation 
effort  is  a further  optimization  of  the  program  so  it  uses  less  computer 

{ , 

time.  The  bulk  of  the  program  execution  time  is  spent  in  updating  the 
wave  modes  so  the  best  way  to  save  time  is  to  reduce  the  number  of  v/ave 
modes  followed  during  an  evolution  calculation  and  to  reduce  the 
frequency  of  mode  update.  The  evolution  results  indicate  the  feasibility 
of  these  reductions  because  only  the  few  fastest  growing  modes  really 
influence  the  ion  distributions  and  because  the  evolutions  can  tolerate 

a small  degree  of  breakdown  without  much  loss  In  the  accuracy  of  the 

I 1 

transport  coefficients. 

There  seems  to  be  no  convenient  way  to  extend  the  program  to 
include  electromagnetic  effects.  It  would  have  to  be  modified  extensive- 
ly unless,  somehow,  the  ions  could  be  assumed  to  be  electrostatic.  This 
extension  would  allow  the  evolution  of  magnetized  distributions  with 
very  high  relative  velocities.  The  electrostatic  criterion,  Eqn  (78), 
is  violated  at  these  high  velocities  because,  in  accordance  with  Eqn  (18), 
the  wave  k becomes  too  small.  The  desirability  of  such  an  extension  is 
not  clear,  however,  because  the  ease  of  the  electron  motions  along  the 
magnetic  field  lines  would  probably  preclude  any  unstable  modes  which 
propagate  at  large  angles  to  the  streaming  direction.  Therefore,  one- 
. dimensional  hydrodynamic  theory,  which  already  has  been  done,  would 

suffice. 
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To  conclude  this  chapter,  a quick  summary  of  the  entire  dissertation 
Is  as  follows:  In  Chapter  I,  the  connection  was  made  between  the 
dissertation  time  evolutions  of  plasmas  with  three  particle  species  and 
the  debris-air  coupling  problem.  In  Chapter  II,  the  quaslllnear  equa- 
tions were  presented,  and  some  physical  discussions  were  made  regarding 
these  equations.  In  Chapter  III,  the  solution  procedure  for  the  equa- 
tions was  explained.  The  solution  of  the  equations  describes  the  slow 
time  evolution  of  the  Ion  distributions  which  Is  caused  by  a spectrum  of 
weak  electrostatic  waves.  In  Chapter  IV,  the  results  of  the  evolutions 
were  presented.  And  finally.  In  this  chapter,  the  results  were  sum- 
marized. The  general  conclusion  Is  that  a fair  amount  of  momentum 
transfer  between  the  debris  and  the  air  occurs  for  the  low  relative 
velocity  plasmas  but  little  occurs  for  the  high  velocity  ones. 

As  a final  note,  recall  that  the  dissertation  program  can  be 
applied  to  other  areas  than  the  high  altitude  nuclear  burst  problem. 

Such  other  areas  Include  the  bow  shock  of  the  terrestrial  magnetosphere, 
the  magnetotail,  and  the  controlled  thermonuclear  reaction  problem.  In 
fact,  the  dissertation  techniques  are  often  applicable  for  a general 
physical  problem  which  Includes  two  low  density  warm  ion  beams  and  a 
single  electron  distribution. 
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Appendix  A 


The  Derivation  of  the  Magnetlc-Fleld-Free 


Juaslllnear  Equations 


The  derivation  of  the  field- free  quasll Inear  equations  for  electro- 


static waves  starts  with  these  nonlinear  Vlasov-ttaxwell  equations: 


y*  + t . y*  = 

vs*  c)  sv 
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ando<  denotes  plasma  species  (debris,  air,  or  electrons),  q^e  denotes  the 
charge  of  the  species,  and  m*  denotes  mass.  The  Initial  density  Is 
assumed  to  be  uniform.  The  Vlasov  equation  is  spatially  averaged  with 


the  average  defined  as 
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Note  the  quantity  Ef^^  Is  of  second  order.  Hence,  the  time  rate  of 
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change  of  Is  of  that  order.  When  Eqn  (A-9)  is  subtracted  from  Eqn 
(A- 1 ) , 
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With  the  spatial  Fourier  transforms  defined  as 
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the  transforms  of  Eqns  (A-9)  and  (A-10)  are 
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The  Green's  function  for  the  operator  on  the  left  side  of  Eqn  (A- 14)  is 

- -JL  (A-15) 

so  the  integral  equation  equivalent  to  Eqn  (A-14)  Is 
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+ dsil #'  i a -*;>')  ■ ^ ^ 


16) 
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When  Eqn  (A-16)  Is  substituted  into  itself  to  iterate,  the  first  few 
terms  of  the  result  are 

* ^ 


^ [/"  ^»-*\>  J -o') 


-J.  A - Ar  Jf  ’ 


] 


(A-17) 


The  first  term  on  the  right  hand  side  of  Eqn  (A-17)  Is  the  free-streaming 
term  based  on  the  Initial  disturbance  to  the  distribution.  The  last 
term  describes  the  effect  of  coupling  of  energy  from  one  wave  mode  at  a 
given  k into  another.  This  term  is  ignored  as  each  wave  mode  is  assumed 
to  proceed  Independently  from  all  of  the  rest.  This  assumption  Is  good 
when  the  wave  amplitudes  are  small,  first  order  quantities.  The 
meaning  of  "small"  Is  described  In  the  Diffusion  Analogy  section  of 
Chapter  II. 

The  electric  field  is  assumed  to  be  In  the  form: 

— -<•  --tf  ( A 


E(A,J-)  - E(A,^~o') 
where  £ ( X y>)  - F 


because  ~ x Iz  (A,* 0 = o 
The  variable  s(k,t)  Is  defined  so 

* CoU,jf) 

c*  - 

or  (X,  = ) etf' 


(A- 18) 
(A- 19) 
(A- 20) 


(A-21) 
(A- 22) 
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The  time  variation  of  to  Is  assumed  to  be  slow  compared  with  where 

co  = to  + i to..  In  fact,  when  to  /i  co( t),  Eqn  (A- 18)  reduces  to 
r i 

^ Co  (X  ) 

(A-23) 

From  the  condition  that  E(j?,t)  must  be  real  so  E*(x,t)  = E(x,t)  and  from 
the  definition  of  the  Fourier  transform, 

E*(-X,Ss>  - (A-24) 

When  Eqn  (A-18)  Is  solved  for  s(k,t)  and  Eqn  (A-24)  Is  used, 

(-XjS}  - (A-25) 


and 


- — co 


a,  a 


(A-26) 


The  Inversion  property  of  4*>{k,t)  In  Eqn  (A-26)  is  often  very  useful. 
The  second  term  on  the  right  hand  side  of  Eqn  (A- 17)  Is  changed  by 

“s,nj  r 

iZ-zA  .*  _ — — 


^~4z.^  - o>(Xj*  or £(I'  ^o)  <a-27> 


and  by  Integrating  by  parts.  The  integrated  part  is  saved,  and  the 

integral  is  neglected  because  it  Is  of  third  order.  Eqn  (A- 17)  then 

becomes  after  the  mode  coupling  term  Is  dropped: 

_»  ^ ^ -J. 

-ut  (fl.28) 

'K1»t  ^ 

Equation  (A-2),  Fourier  transformed,  is 

(A-29) 


100 


AFIT/DS/PH/78-1 


where  Eqn  (A-20)  is  used  again.  When  Eqn  (A-28)  Is  inserted  into  Eqn 
(A-29), 


/ 


- I 


•nr-ti  S 


r 

( J -*  ^ 

> r -^rx,s)] 


= J 
«(. 


-t  -e. 


-t_-A  £(J£j  X - <0 


J ^~°X. 


^ *a^l±L  r /-> 

s***  Ji?  X r T - 


(A- 30) 


where,  as  necessary,  the  integrals  are  taken  below  the  singularity  in 
the  denominators  to  Insure  the  results  are  analytic  across  the  real  axis 
in  the  complex  frequency  plane.  The  Integration  symbol  J denotes  this 
process  because  the  process  Is  equivalent  to  integrating  along  the  Landau 
contour.  The  first  and  second  terms  on  the  right  hand  side  of  Eqn  (A- 30) 
are  the  free-streaming  terms.  These  terms  are  not  associated  with 


growing  or  decaying  waves  and  phase  mix  to  zero  in  times  of  the  of 

-l 


As*r) 


(A- 31) 


where  z^v  is  the  width  of  ^j(k,v,t=0)  for  the  first  term  and  is  the 
width  of  £0(v,t)  for  the  second  term.  These  are  the  times  necessary 
for  the  particle  free  streaming  to  smooth  out  Initial  spatial  inhomo- 
geneities. The  time  asymptotic  form  of  Eqn  (A-30)  is,  therefore, 

^ r j,. . t ■ -= 


i-i 


a. 


•<  u 


fJ 


Jbv-  — =3- 


— O 


(A- 32) 


[_Ja.>  AT  ~ 

which  Is  the  linear  dispersion  relation  for  = Co( k,t).  This  equation 
Is  valid  for  both  damped  and  growing  wave  modes.  Mote  however,  the  waves 
cannot  be  too  strongly  damped  because  then  the  waves  will  die  out  before 
the  free-streaming  terms  go  away.  The  damping  time  which  Is  approximate- 
ly^^ must  be  greater  than  the  time  given  by  Eqn  (A-31). 


AFIT/DS/PH/78-1 


The  quasi linear  equation  for  the  time  rate  of  change  of  the  distri- 
bution results  when  Eqn  (A-28)  Is  Inserted  Into  Eqn  (A-13): 

) _ __  # 


~J.Jl  ‘s^r-sf 


— oO 


- f***-  E (X t ■£'}  — -&- £ (yhjjt  - *Q  (A— 33) 

y<  [ M -yy  — jZyt)  ] e>y^r  J J) 


/Hi  . 


When  the  free-streaming  terms  are  neglected  and  the  integral  Is  taken 
below  the  singularity, 

('V'jJt)  _ #«*  x ^ f //  f J_  — l — Ff-X  /\  p ( £ y\ 


y.  y 

X*-  ^ 


(A- 34) 


The  free-streaming  terms  can  be  neglected  for  times  greater  than  Ukv) 
where  ^k  Is  the  width  of  the  wave  spectrum  and  visa  typical  particle 
velocity.  This  criterion  is  similar  to  Eqn  (A- 31) . 

#o  -A 

The  spectral  wave  energy,  £ (k,t),  and  the  quasll inear  diffusion 

coefficient  for  species*,  D^v.t),  are  respectively  defined  as 

e / T = jUh.  _ (. _ £(~^ij  £(£>#) 

fc  ^ y r-Lfr)1  srft~ 


-1 


(ergs)  (A-35) 


x x 


f*i  -*L  r .-f  £ (■£■,! ^ 


(A- 36) 


Therefore,  Eqn  (A-34)  becomes 


, A • [v^  >v  %a£p£>]  tft-37) 

dAr  L ‘ 4*  \ 

The  time  dependence  of  the  spectral  wave  energy  from  Eqns  (A-18),  (A-24), 
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7,  X •**  £ -4  ( ^ 


(A-25),  and  (A- 35)  Is 

When  Eqn  (A- 38)  Is  differentiated  and  Eqn  (A-21)  Is  used. 


(A- 38) 


3 c 

O'* 


(A- 39) 


A nx>re  convenient  form  of  the  diffusion  coefficient  results  when  Eqn 
(A-36)  Is  broken  Into  its  real  and  Imaginary  parts  and  Eqn  (A-26)  is 
used: 

- fir  ^5  J (A',0) 


where  £ (~^j  jf  - o')  = £ - o')  (A-41) 

has  been  assumed. 

Equations  (A-32),  (A-37),  (A-39),  and  (A-40)  are  the  simultaneous 
set  of  quasillnear  equations  for  wave-particle  interactions. 
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Appendix  B 

Conservation  In  Quaslllnear  Theory 


The  equations  of  quasllinear  theory  conserve  natter,  momentum,  and 
energy  as  can  be  seen  by  taking  the  first  three  velocity  moments  of  Eqn 
(10).  The  first  moment  is 


(B-l) 

(B-2) 


- o 


(B-3) 


The  surface  Integral  Is  at  Infinity  because  the  velocity  Integral  Is  over 
all  space.  Equation  (B-3)  shows  that  the  species  density  Is  a constant 
so  particles  are  conserved.  The  next  moment  is 


where  the  chain  rule  and  the  divergence  theorem  have  been  used.  When 
Eqns  (9)  and  (A- 36)  and  the  sum  of  Eqn  (B-5)  over  all  of  the  species  are 
all  combined. 


= o 


(B-6) 

(B-7) 


Equation  (B-7)  follows  because  the  Integrand  Is  odd.  Therefore,  momentum 
Is  conserved.  The  third  moment  Is 
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— e>6  — 


(B-8) 


J>_ 


sptc'ics  U 
pari'  ct< 
<x<rqy 

val 


- . 4-  f ^ f 

-<-  ^>v\^  / V - CO^Jx 


(B-9) 


where  as  before  Eqn  (A-36)  Is  used.  Equation  (B-9)  is  broken  into  two 


separate  integrals  one  of  which  Is  odd  by  letting 

» -4 

^ ^ — co  -t  to  (B-10) 


The  symmetry  property  of  co,  Eqn  (A-26),  and  Eqn  (9)  are  also  used. 
Therefore,  after  the  sum  of  Eqn  (B-9)  over  all  of  the  species  is  taken. 


Aj 

d* 

‘ f.-i.f  , * 

par-ii c l< 

<7*» 

~ ~ ^ j ^ 

- 0* 

(B-ll) 

i v.l  . 

r - J dA 

(B-12) 

I*  (va«  J 

= - J-  \ iurjz  / 

[_  v«!  J 

(B-13) 

A loss  in  particle  energy  is  balanced  by  a gain 

in  wave  energy  so  total 

energy  is  conserved. 
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Appendix  C 

Electron  Dielectric  for  a Nonzero  Magnetic  Field 


The  dielectric  is  derived  from  the  first  order  perturbation  of  the 

electron  distribution  function  calculated  in  the  Lagrangian  coordinate 

' \ 

system.  The  Lagrangian  system  follows  the  zero  order  trajectories  of 
the  electrons.  The  derivation  is  very  long  and  tedious  so  only  the 
highlights  are  presented.  For  further  details,  see  Ref  16. 

The  rate  of  change  of  the  distribution  along  the  trajectory  "x  = 


x(t)  Is 

The  change  along  the  zero  order  trajectory  Is  given  by 


f _ — 


* — t ", 


F&l  = zero  order  solution  = 0 

- ^eo  (A?)  t- 

e = lei 


(C-l) 


(C-2) 

(C-3) 

(C-4) 

(C-5) 

(C-6) 


Equation  (C-4)  describes  the  helical  motion  of  a charged  particle  in  a 
uniform  magnetic  field.  Vlasov's  equation,  Eqn  (1),  after  linearization 
Is 


The  result  when  Eqns  (C-2)  and  (C-7)  are  compared  is 
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rate  of  change  of  ~|  _# 

( - f , along  the  zero / - ±.  ( £ + ^ *!•).  ^At  (C-8) 

l It  / [order  trajectory  J c ' ** 

where  Eqns  (C-4)  and  (C-5)  have  been  used  to  eliminate  (*r  * *0* 

The  solution  of  Eqn  (C-8)  for  growing  wave  modes  is  (The  solution 


for  damped  modes  is  done  by  analytic  continuation.) 

Jt 


^ f S')  + — 

— c+ 

where  the  integration  is  done  along  the  zero  order  trajectory.  Equation 
(C-9)  is  rewritten  using  Maxwell’s  equations  assuming  plane  wave 


sol utions: 


/«»  (*r,  t')  - ^ ^ctf' 


U.  A*-  • s* 


r • x 'A  ^\r 

L * ~Z 


Z — J‘  (c-,o> 


This  equation  in  the  Eulerian  (fixed)  coordinates  v\  t,  and  t is 


^ - a - 


•t  + 


(-  SLjt  'T'  + M ( I - c-o-a  Sl+T  + 


JiAT 


— sbsx  ( I C<ra 


SL<T'\)  + ^ ^ 


^ f T"  "f  **^*jft 


' [/“>-<-  * (/*<>*  /«i)]  + £2  [ ~Zr(--^*r  r 
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+ •^Ki  ^ -1*  •****  -H-«. ^ 


where 


_SIA=  electron  gyro  frequency 
—»■  A 

Ba  z Sc  * 


uiiaj 


/-  JL  d/eo  - L ^/*° 

*0±  ' ^ d/V*  ~ sv^ 

/«*  = ^ 


(C-ll) 


(C-12) 


(C-13) 


(C-14) 


{C— 15) 


(C- 16 ) 


The  last  two  equations  follow  from  the  azimuthally  symmetric  form  of  the 


zero  order  distribution; 


syri  x 
~ ~f*»  ( -3>  3 y 


(C— 17) 


where  yyr^  - sy/£  + ^ 

Note  Equation  (C-ll)  Is  of  the  form: 


- -c- 


(.SVj  &*0 


(C-18) 


(C-19) 


The  remaining  equations  will  only  deal  with  the  Fourier  amplitude, 
fei(v,k,c°).  The  phase  factor  Is  dropped. 

The  spatial  average  of  fej(v,k,t*>)  in  the  direction  perpendicular 
to  the  ambient  magnetic  field  Is  defined  as 


e*-» 


(C-20) 


The  zero  order  distribution  Is  assumed  to  be  a Maxwellian  In  the  perpen- 


AFIT/DS/PH/78-1 


dicular  direction  so 


a.  n.'iry:TJtl 


(C-21) 


where  the  distribution  In  the  parallel  direction  has  yet  to  be  specified. 
The  spatial  average  of  Eqn  (C-ll)  In  terms  of  rotating  coordinates  Is 


/ / N . -s-x-aZ  \*r 

<T->i  '*"■*  ^ = . 


JlA  -f  co  — -»A  ) "T" 


where 


~ ±.  Lx  '2- 

* ‘ a -Ai J>l* 

Lt = average  gyro  radius 


* ^y- 


h*  S~  + -£ 


(C-23) 


(C-24) 


(C-25) 


(C-26) 


W 


9-«o  C^y 


l CO 


40 

T + XT, 


<i  y 


-^**y  <3. 


Xh-  x.n)  * 


Sa  = ^ 


(C-27) 


(C-28) 


(C-29) 


The  Bessel  functions,  On (1 ^ ) and  In(  ^),  In  Eqn  (C-28)  are  respective- 
ly the  ordinary  and  modified  Bessel  functions  of  the  first  kind.  After 
the  distribution  In  the  z direction  Is  specified  as  the  following 


Maxwellian  with  a zero  drift  velocity: 

r ...  r//»! 


yV~ 2 


-[tffxTi#]  *X7U// 


(C-30) 
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r 


the  total  spatial  average  of  the  first  order  distribution  becomes 


— 

_ -A  -A  n 

o<» 

•2>j 

y>»*  — «o 

-fo 

(I:-')'  £ 

] F,  6>u) 

i7a 


tia. 
( j~. 


where 


•V 


SCO  « - - %(*„') 

^)  = /^  + «u.  r-M] 


Co 


+ /vi  j 

~AX  M **7i 


// 


*2  (*u)  9 plasma  dispersion  / r J* 

* J function  (Ref  17)  - — ' ^ 


r ^ 

i*- 


frF  ). 

The  electric  field  Is  assumed  to  be  electrostatic  so 
E cJ)  = J*  <f>  (Jt-)  co) 


(C-31) 
(C-32) 
(C- 33) 
(C-34) 
(C-35) 

(C- 36) 


(C— 37) 


where  4 Is  the  electrostatic  potential.  When  Eqns  (C- 37)  and  (C-31) 
are  combined  with  the  Fourier  transformed  Poisson’s  equation, 

o * 

E (Xj  *0  = X ( C- 38 ) 


the  following  form  of  the  electron  dielectric  results: 


6.  - 


Co 


, c,. 


VO  - T7 


■|  tSc,,  A [ ( X.e  (*,) 
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+ (co  --SU)  f 

(C-39) 

where  °°f>*  c P^asma  frequency  = "V 

(C-40) 

(C-41) 

Zi  r Z// 

(C-42) 

/x  T 

= dermal  speed  - J *'tj  * 

' f 

(C-43) 

The  Infinite  sum  has  been  eliminated  by  assuming 

X i 

(C-44) 

so  Te  O')  - 1 

(C-45) 

= -kr 

(C-46) 

TM2x(X)  Is  of  order  X or  higher 

(C-47) 

Only  the  zero  and  first  order  terms  In  X are  retained  In  Eqn  (C-39). 
This  result  Is  also  valid  for  damped  modes  because  the  definition  of  the 
plasma  dispersion  function  includes  the  analytic  continuation. 

Equation  (C-39)  appears  to  be  singular  for  propagation  exactly 
perpendicular  to  the  ambient  magnetic  field.  However,  the  lim.t  of  the 
dielectric  as  kz~ *0  Is 


which  Is  well  behaved. 


Sli 


(C-48) 
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Appendix  D 

Analytical  Form  of  the  FI  eld- Free  Dielectrics 


An  analytical  form  of  the  field-free  dielectric  Integrals  In  Eqn 
(9)  results  when  the  following  two-temperature  Maxwellian  Is  used: 


/ / f I l"'1  (d.: 

where  the  streaming  velocity,  V^,  Is  assumed  to  be  In  the  x direction. 
To  begin  the  derivation,  Equation  (9)  Is  rewritten  as 


(D-l) 


\ 2 
where  co^  = (species  plasma  frequency)  * 


'*o 


(D-3) 


— I 

-'viw# 


(D-4) 


When  Equation  (D-l)  Is  Inserted  Into  Eqn  (D-4)  and  the  result  Is 


multiplied  by  k. 


'X'Tr 


\z  — 


(D-5) 


where  the  one-dimensional  root-mean-square  thermal  velocities  In  the  x 


and  z directions  are  defined  as 


C*'X,  2 


2 


(D-6) 


The  most  convenient  way  to  do  the  integral  1^  In  Eqn  (D-2)  Is  to 
perform  a coordinate  rotation  In  the  v -v_  plane  to  align  the  new  v 
axis  (vxi)  with  the  wave  propagation  direction.  Then,  the  singularly 
In  the  denominator  appears  only  in  the  v , coordinate  so  the  integration 
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In  v2,  can  be  done  directly.  The  relationships  between  the  wave  k and 


velocity  vectors  In  the  old  and  new  systems  are 


(;:) 


A* 

> u. 


6'0-si  & 


'Ur 


(D-7) 


When  these  relations  are  used  along  with  the  fact  that  the  Jacobian  for 
a rotation  Is  unity,  the  result  after  much  algebra  Is 

/JTI  i f,  , fc^^'  * CcJy -<^0  1 

-t*  ' 'I'Tr  tJtsyCjlljJ*'*'**'  1 — Vgl 


to  - 


p*  m 

f <v\^'C4-a6>  T -cjj)*] 


( * c,\  c4 

where  c**x&  ■+  cj\^sx9 

~ f c*n.'L& 


(D-8) 


(D-9) 


(D-10) 


The  angle  6 Is  that  between  k and  the  vx  axis.  After  the  integral  In 
v2,  Is  done  and  the  result  Is  simplified,  the  following  simple  expression 


appears 


l C / (sVy'  - V-tC+iB) 


1"  r 13t  cje  (A 


V — . . a 7T  -B- 


(sV^c'  - ^ 

* c?& 


(D-ll) 


" w ^ - « W[Xi 

This  form  Is  related  to  the  derivative  of  the  plasma  dispersion  function 
as  defined  by  Fried  and  Conte  (Ref  17).  This  derivative  Is 


tD-12) 

Hence,  the  final  form  of  the  dispersion  relation  results  when  Eqns  (D-2), 


(0-11),  and  (0-12)  are  combined: 


C*JmA 


(D-13) 
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where  Is  the  root-mean-square  thermal  width  In  the  © direction  from 
Eqn  (0-9)  and 

WiXl  - 


nJx  c 


(D- 14) 


Note  coi7c^ecan  be  rewritten  in  terms  of  the  species  Debye  length  In 


P* 

the  & direction  because 

* 


(0-15) 


Incidentally,  the  dispersion  relation,  Eqn  (0-13),  has  the  same  form  as 
a one-dimensional  dispersion  relation. 

The  dispersion  relation  is  valid  for  all  values  of  the  complex 
frequency  because  the  definition  of  the  plasma  dispersion  function  con- 
tains the  analytic  continuation  across  the  cor  axis  Into  the  lower  half 
of  the  complex  frequency  plane.  This  continuation  Is  equivalent  to 
conducting  the  Integration  in  Eqn  (0-12)  along  the  Landau  contour. 
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I 

'i 

Appendix  E 

The  Derivation  of  the  Anomalous  Transport  Coefficients 

The  calculation  of  the  anomalous  transport  coefficients  for  momen- 
tum and  energy  transfer  Is  done  using  the  formulas  derived  below.  These 
formulas  are  obtained  from  the  first  three  velocity  moments  of  the 
spatially  averaged  Vlasov  equation,  Eqn  (A-9).  This  process  gives  the 
standard  fluid  equations  of  continuity,  momentum  transfer,  and  energy 
balance  with  anomalous  collision  terms  Instead  of  the  normal  collision 
terms  due  to  binary  interactions. 

In  the  derivation,  the  following  definitions  are  used: 


^*0)  * Jjj?  jj. 

(E-l) 

'*'*t(‘*/*}  9 X/1*  *4**  I*/'*' j*) 

(E-2) 

**- 

(E-3) 

o° 

’/n  V)M  (*,  *)  = ^ 

(E-4) 

= Jotc  (^rJ 

(E-5) 

where  T^  Is  in  degrees  Kelvin,  and 

= random  velocity  = (E-6) 

For  convenience,  Eqn  (A-9)  is  repeated  here: 

/£.  fVW  \ (E"7) 

' s***  > oSr  / 

The  velocity-to-the-zero-power  moment  simply  gives  the  equation  of 
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continuity  which  is  not  of  interest  because  the  plasma  is  assumed  to  be 
initially  homogeneous.  The  next  two  moments  are  done  using  the  methods 
in  Seshadri  (Ref  7:  27-45). 

The  velocity- to- the- first- power  moment  of  Eqn  (E-7)  is 

33.J-V.-~v7,]  = - ^ (E-8) 

= (E-9) 

where  the  chain  rule,  the  divergence  theorem,  and  Eqns  (E-2)  and  (E-3) 
have  been  used.  Equation  (E-9)  is  the  momentum  transfer  equation  with 
the  right  hand  side  representing  the  anomalous  change  in  momentum  density 
due  to  wave-particle  interactions. 

The  third  moment  of  Eqn  (E-7)  Is 

rt  •>>,  (*-Y*?]y  <E-»> 

where  the  above  procedures  and  Eqn  (E-5)  have  been  used.  Continuing, 

<b'C  (i  ~ (E-13) 

where  Eqns  (E-2)  and  (E-4)  have  been  used.  Equation  (E-13)  is  the 
energy  or  heating  equation  with  the  right  hand  side  giving  the  anomalous 
heating  rate. 

Equations  (E-9)  and  (E-13)  are  put  into  a more  convenient  form  by 
using  the  Fourier  transformed  Poisson's  equation,  Eqn  (A-29), 
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iA  £ (A,  Z)  - WlL  pi-*- 


(E-14) 


and  the  dispersion  relation. 


/ = Z«. 


(E- 15) 


The  plasma  dielectric,  6^  , from  Eqn  (9)  is 
..sn-XI  "7* 

*1%,  u , . x-  -t-izr1- 


€ - \<&— ^rr 

* - «o(Z,A] 


(E-16) 


When  both  sides  of  Eqn  (E- 15)  are  multiplied  by  ikE(lc,t),  Eqn  (E-14) 


used,  and  corresponding  terms  of  the  sums  equated; 


"'ll  A (£> 'Zj*}  ~ Hlt'f+JL 

This  equation  along  with  the  definition  of  spectral  wave  energy,  Eqn 
(A-35) , are  the  main  equations  used  to  recast  the  anomalous  collision 
terms.  In  addition  to  these  equations,  the  Integral  form  of  the  delta 


(E-17) 


function. 


£ C A A ^ “ (^y  j <Ax  -4- 


u>z\ 


(E-18) 


and  the  symmetry  of  the  dielectric. 


(E-19) 


are  used.  Equation  (E-19)  follows  from  Eqn  (E-16)  broken  into  its  real 
and  imaginary  parts  and  from  Eqn  (A-26).  Therefore,  from  Eqn  (E-9), 


— 

■ Sfj’  s 


(E— 20) 
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= $dl.X  €(XtX)  £*  (u>  (X,/'^/  X^  (E-2i) 

= 2 j tU\  Jl  £(X,X)  ( C^XJ^)JX)J  (E-22) 

and  from  Eqn  (E- 13) , 

~T^<»  ‘ (s* (E-23) 


= *,<<£>■  v xu>o>:>^,  j & ^ iftsp 

-*• 


~ dJi  £(■£.)"$}  oo  to 

— » 

o* 

- 2 Z'V^€tli(l~(Zj'>j£)  (E-25) 

= -n^dA  £ (2,/)$c*^(X,/)  fal^CcufX,/);  X)] 

— c*>  ^ 

+ (E— 26) 

where  In  the  first  term  on  the  right  hand  side  of  Eqn  (E-23),  Eqn  (A-28) 
without  the  free  streaming  terms  and  Eqn  (B- 10)  were  used  to  obtain  Eqn 
(E-26).  Equations  (E-22)  and  (E-26)  are  the  normal  fluid  momentum 
transfer  and  energy  balance  equations  for  a collislonless  homogeneous 
plasma  with  added  anomalous  collision  terms  on  the  right  hand  sides. 
Because  the  plasma  Is  homogeneous,  all  of  the  changes  In  the  momentum 
and  energy  densities  are  due  to  anomalous  transport. 

Although  the  above  derivation  assumes  the  plasma  Is  magnetic  field 


118 


AFIT/DS/PH/78-1 


free,  the  generalization  to  a plasma  with  a uniform  magnetic  field  is 
very  straight  forward.  In  fact,  the  resulting  anomalous  transport 
coefficients  are  the  same  as  those  in  Eqns  (E-22)  and  (E-26).  With  the 
addition  of  the  magnetic  field,  Equation  (E-22)  becomes 

— a4 

where  Bq  is  the  magnetic  field.  The  energy  balance  equation  is 
unchanged. 
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Appendix  F 

Numerical  Form  of  the  Quasllinear  Equations 


In  this  appendix,  the  quasllinear  equations  are  each  cast  Into 
their  numerical  forms  by  the  extensive  use  of  symmetries.  The  resulting 
equations  are  the  ones  put  Into  the  solution  program  controlled  by 
I0HI0N.  The  diffusion  equation  is  discussed  first. 

The  diffusion  equation  is  solved  only  on  the  upper  half  of  the 
computer  solution  velocity  plane  shown  In  Fig.  6 because  the  equation  Is 
symmetric  with  respect  to  reflection  across  the  vx  axis.  This  half 
plane  is  mapped  onto  a rectangle  (the  computer  mesh)  with  Its  lower  left 
and  upper  right  corners  at  (-1,0)  and  (1,1)  respectively.  This  mapping 
is  described  by  the  following  relations: 


where  - I - - I 


o £ ^£/i3  £ J 


SV'sro  ***** 

’ — 1/ 

Vt>  A 

(F-l) 

-Jl  - 

(F-2) 

(F-3) 

(F-4) 

The  variables  u^  and  uz  are  the  mesh  coordinates.  Three  separate  sets 
of  the  three  transformation  parameters;  vXQ,  vZQ,  and  Vcft;  are  used  for 
each  of  the  plasma  species;  debris,  air,  and  electrons.  Therefore,  a 
given  value  of  a mesh  coordinate  refers  to  three  separate  velocities 
depending  on  which  species  is  being  referenced. 

The  transformation  parameters  are  determined  from  the  initial 
species  Maxwell ians  and  are  held  fixed  throughout  the  time  evolution. 
The  parameter  Is  chosen  so  the  computer  mesh  Is  centered  on  the 
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electron  distribution,  or  for  the  ions,  is  centered  one  thermal  width 

in  the  vx  direction  away  from  the  distribution  peak  on  the  side  toward 

the  other  distribution.  The  scale  parameters  v and  v are  respec- 

xo  zo 

tively  set  equal  to  three  times  the  species  thermal  width  in  the  vx 
direction  and  that  in  the  vz  direction.  These  settings  of  the  trans- 
formation parameters  cause  the  majority  of  the  mesh  points  to  correspond 
to  the  velocities  of  the  interesting  parts  of  the  species  distributions. 

For  good  numerical  results,  the  mesh  spacing  (the  number  of  mesh 
points  used)  in  the  ux  and  u2  directions  are  chosen  so  at  least  a few 
mesh  points  correspond  to  the  resonant  region  (where  Eqn  (18)  holds). 

The  mesh  spacing  is  also  chosen  to  Insure  sufficient  accuracy  of  the 
dispersion  relation  integrals  described  later.  This  spacing  is  fairly 
critical  because  the  mode  growth  rates  depend  on  the  distribution 
slopes. 

An  interesting  result  of  the  above  mesh  transformation  procedures 
is  that  a given  plasma  species  appears  the  same  on  the  mesh  regardless 
of  its  temperature,  mean  velocity,  or  density.  For  example,  a sketch 
of  the  initial  debris  Maxwellian  is  shown  in  Fig.  F-l.  Therefore,  the 
numerical  accuracy  of  the  evolution  is  relatively  unaffected  by  changes 
in  the  plasma  distribution. 

The  diffusion  equation  for  growing  modes,  after  expansion  of  the 

dot  products  and  conversion  into  the  mesh  coordinate  system,  is 

.3 


d*  = — r— 

+ 'x  /a  (^,*2 

■ IN  f 
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COMPUTER 

MESH 


DEBRIS  THERMAL  WIDTH 


Fig.  F-l.  Initial  Debris  Maxwellian  Distribution  on  the 
Computer  Mesh. 

• ^ ></t»  ^ * ~ ^ P C ✓ _ .X 
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^io 


.ilZ^lL  - D (M  f)  l^i(  (F-5) 

where 

^ - ^U^&ttckXrM^  <-> 

. N •v'fT  /7  JlSK 

l~>  ■)0='1  “^[  (X.^U  < F- 9 > 

6««  (~^---,a^j)‘ul(Z.^-~4Z.*tt  (F-10) 
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Equations  (F-9)  and  (F-10)  resulted  from  the  combination  of  the  velocity 
derivatives  of  the  diffusion  coefficients  given  by  Eqns  (F-6)  to  (F-8). 
Equation  (F-5)  appears  to  be  an  elliptical-parabolic  equation  with  two 
spatial  variables,  ux  and  uz.  The  right  hand  side  Is  elliptic  because 
the  coefficients  to  the  second  order  terms  are  positive  definite,  and 


/ 

£ 

* 

c 

r-  -c  11 

L -n.  J 

For  damped  modes  or  a combination  of  damped  and  growing  modes,  Eqn  (F-5) 
must  be  analytically  continued.  This  process  is  numerically  very 
difficult  and  Is  not  necessary  because  the  growing  modes  dominate  the 
solution  for  the  times  of  Interest.  (The  fastest  growing  nodes  build 
up  the  largest  wave  energies  which  modify  the  distribution  the  most.) 

A mixture  of  damped  and  growing  modes  appears  only  late  In  the  calcula- 
tion anyway. 

Because  the  plasma  species  distributions  all  go  to  zero  as  vx  or 
vz  goes  to  infinity,  the  boundary  conditions  for  the  diffusion  equation 
are  relatively  simple: 


^olt,  ~ 1 J jjT)  ~ 0 

(F-12) 

■ft*  ~ ° 

(F-13) 

C j ~ O 

(F-14) 

The  symmetry  of  the  upper  and  lower  half  planes  determines  the  boundary 
condition  on  the  mesh  bottom.  The  initial  condition  for  each  distribu- 
tion Is  the  Initial  Maxwellian, 

(sv^-V*)2- 

( F- 15 ) 
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where  vx  and  vz  are  given  by  Eqns  ( F-l ) and  (F-2). 

In  the  method  of  lines,  the  spatial  derivatives  on  the  right  hand 
side  of  Eqn  (F-5)  are  written  In  terms  of  finite  differences,  and  the 
resulting  set  of  simultaneous  ordinary  differential  equations  Is  evolved 
In  time  by  the  GEARB  package.  The  interior  Including  the  bottom  of  the 
mesh  Is  broken  up  Into  a series  of  vertical  lines  of  mesh  points  with 
the  points  sequentially  numbered.  The  resulting  numbering  scheme  Is 
shown  In  Fig.  F-2  for  a nx  by  nz  mesh  where  nx  and  nz  are  respectively 


the  numbers  of  mesh  points  in  the  ux  and  uz  directions.  A finite 
difference  equation  Is  written  for  each  of  the  numbered  points  using 
the  following  finite  difference  formulas: 


" 

c)^Ao 

(F-16) 

ft:  1 7 \ 

V -A.,*. 

\F-17 ) 

t 

*1 

^>1 

-j Oft  -7^  -t 

( F— 18 ) 

At 

(F-19) 

T - . / 

7, -A ./Y 

1 -At  '/ 

(F-20) 

7 X* 

where  hx  and  hz  are  the  uniform  mesh  spaclngs,  and  1 Is  the  mesh  point 
of  Interest.  As  shown  in  Fig.  F-2,  this  is  a nine  point  differencing 
scheme.  For  the  difference  terms  on  the  top  and  sides  of  the  mesh,  the 
boundary  conditions,  Eqns  (F-12)  to  (F-14),  are  used.  The  equations  for 
the  bottom  mesh  points  use  the  symmetry  of  the  distribution.  For  these 
points,  Eqns  (F-16)  and  (F-19)  are  unchanged,  Eqns  (F-17)  and  (F-20) 
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vanish,  and  Eqn  (F-18)  becomes 


* | _ i ^ */-* 


AX 


(F-21 


'*  Hotter  ^ 

The  finite  difference  equation  coefficients  are  always  evaluated  at  the 


point  of  Interest. 

The  equation  for  mesh  point  1 from  Eqns  (F-5),  (F-12)  to  (F-14), 
and  (F-16)  to  (F-21)  Is 

3L 
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The  set  of  simultaneous  ordinary  differential  equations  each  given 
by  Eqn  (F-22)  or  (F-23)  Is  the  set  which  Is  solved  by  the  GEARB  package 
using  backward  difference  formulas.  The  set  Is  stiff  because  the 
coefficients  contain  the  exponentially  growing  spectral  energies  of  the 
wave  inodes  and  because  the  method  of  lines  usually  results  In  a stiff 
set  of  equations  (Ref  23).  The  set  Is  of  the  form: 

FJ  = J FJ 


(F-31) 


where 


F J - 


U*  1 
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(F-32) 


F J = 


f A 1 
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(F— 33) 


The  system  matrix  and  Jacobian,  J,  are  the  same  and  are  defined  as 

^ FJ ' 


j-y  - 


3FJ/ 


(F-34) 


The  matrix  and  Jacobian  have  the  banded  structure  shown  In  Fig.  F-3. 

This  banded  structure  Is  exploited  In  the  GEARB  routines  to  save  computer 
storage  and  to  Increase  computational  speed.  The  total  number  of 
equations,  NJ,  Is 
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integrals  are  done  numerically  using  a wave  k grid  determined  during 
the  solutions  of  the  dispersion  relation  for  the  Initial  plasma  distri- 
bution. At  that  time,  the  location  of  the  grid  in  the  first  quadrant  of 
the  k plane  is  determined  by  locating  the  area  of  suitable  unstable 
growing  wave  modes.  The  boundaries  of  the  grid  are  constrained  so  that 

O < J**  < (F-37) 

0 £ -<4*  < °<>  (F-38) 

The  kz  axis  must  be  avoided  because  the  dispersion  relation  used  is 
singular  there.  The  grid  is  chosen  to  be  as  small  as  possible  so  only 
the  significant  wave  modes  are  considered. 

The  first  quadrant  grid  is  sufficient  to  do  all  of  the  integrals. 
The  waves  propagating  in  other  directions,  and  thereby,  having  wave  k's 
in  the  other  quadrants  are  Incorporated  into  the  integrals  by  the  use 
of  these  two  symmetries  of  the  diffusion  coefficient  integrands:  (1) 

The  integrands  are  symmetric  with  respect  to  inversion  in  k,  and  (2) 
the  integrand  in  the  fourth  quadrant  is  that  in  the  first  with  the  signs 
of  the  kz  terms  reversed.  These  symmetries  are  rigorously  true  for 
unmagnetized  ions  and  electrons.  They  also  hold  approximately  when  the 
electrons  are  magnetized.  The  symmetries  result  from  the  physical 
equivalence  of  waves  propagating  in  the  e,  - e , S+IF,  and  - & + TT 
directions.  Thus, 

~ CO'S-AvJ-JLly£)  - - ~ (F-39) 

} { p- 40 ) 

s - E ( F-4i ) 


] 

j 

1 

j 


129 


AFIT/DS/PH/78- 1 


= <£U,A^  (F-42) 


where  the  symmetry  of  the  species  distributions  about  the  vx  axis  and 
Eqn  (A-26)  have  been  used.  These  symmetries  are  also  used  in  the  con- 
version of  the  dispersion  relation  into  its  numerical  form. 

Because  the  frequencies  and  growth  rates  are  assumed  to  be  constant 
during  each  phase  of  the  evolution,  the  only  time  variation  in  the 
diffusion  coefficients  during  a given  phase  is  through  the  spectral 
energy,  £ (k,t).  Equation  (11)  can  then  be  easily  integrated  to  give: 


ec.4,/)  = 


'l  (-A,  4 ) (J  ) 

JL 


(F-43) 


where  tQ  is  the  beginning  time  of  the  current  phase  of  the  evolution 
and  £(£,t=to)  is  the  spectral  energy  accumulated  from  the  previous 
phases.  This  result  is  Inserted  directly  into  the  diffusion  coeffi- 
cient integrands. 

An  initial  thermal  spectral  wave  energy,  £(k,t=0),  is  used  to 
begin  an  evolution.  This  energy  is  assumed  to  be  constant  over  the  k 
grid  and  equal  to  the  plasma  thermal  energy  divided  by  the  available 
volume  in  space  for  electric  fields,  the  Debye  sphere.  (A  point  charge 
is  unshielded  by  the  rest  of  the  plasma  inside  the  Debye  sphere,  and 
hence,  random  electric  fields  can  exist  only  within  that  volume.) 
Therefore, 


-i-xr. 


[o(  Kia) T 


OcA 


m-2. 7 


Od  **wnr 


( F— 44) 


I 


I 


■Iz2r)  (F-45) 


where  the  ak's  are  the  lengths  of  the  grid  sides  and  the  factor  of  one 
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quarter  accounts  for  the  energy  In  the  other  quadrants.  The  exact  value 
of  this  small  thermal  wave  energy  Is  not  important.  In  fact,  it  Is 
usually  increased  a few  orders  of  magnitude  in  the  beginning  of  the 
evolution  calculation  to  save  computer  time. 

Two  separate  procedures  are  used  to  evaluate  the  dispersion 
relation  dielectrics.  The  first  is  used  only  on  the  Initial  Maxwellian 
species  distributions  and  Is  the  faster  of  the  two.  In  this  procedure, 
the  analytical  form  of  the  plasma  dielectrics  given  by  Eqn  (20)  or  (32) 
Is  used.  The  plasma  dispersion  function  Is  numerically  computed.  In 
the  second  procedure,  the  Integrals  in  Eqn  (9)  are  numerically  evaluated 
for  evolved  non  Maxwellian  ion  distributions.  Note  the  analytical 
procedure  can  always  be  used  for  the  electron  dielectrics  because  the 
electrons  are  adiabatic. 

The  numerical  Integrals  In  the  ion  dielectrics  are  done  assuming 
that  co.j (k,t)  f 0.  The  Integral  for  species*,  1^,  Is 
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The  angle  between  k and  the  vx  axis  is  The  last  term  of  Eqn  (F-46) 

Is  the  analytic  continuation  for  negative  growth  rates.  It  is  derived 
using  the  residue  of  the  Integrand  numerator  In  Eqn  (9)  at  the  singular- 
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ity  assuming 
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The  resulting  residue  Integral,  1^,  Is 

~£o<Jl  ~ <&  Ji  • (F~49) 
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It  reduces  to  the  form  In  Eqn  (F-46)  after  the  coordinate  rotation 
procedure  in  Appendix  D Is  used.  This  rotation  allows  the  delta 
function  to  be  Integrated  In  the  vx  direction  of  the  rotated  frame  of 
reference. 

Equation  (F-46)  after  the  coordinate  transforms,  Eqns  (F-l)  and 
(F-2)  are  applied  Is 
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The  first  three  terms  on  the  right  hand  side  of  Eqn  (F-50)  correspond 
to  the  integrals  over  the  upper  half  of  the  velocity  plane,  and  the 
second  three  correspond  to  those  over  the  lower  half.  The  signs  In 
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Eqn  (F-51)  refer  respectively  to  the  upper  and  lower  half  planes. 
The  numerical  version  of  the  Ion  dielectrics  Is,  therefore, 


^ = 


H-'tr- 
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X.  'i. 


(F-52) 


with  1^  given  by  Eqn  (F-50).  With  this  equation,  the  conversion  of  the 
quaslllnear  equations  into  forms  suitable  for  numerical  solution  Is 
complete. 
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Appendix  G 
The  Solution  Program 

The  quasilinear  equations  for  the  ion-ion  instability  in  the  form 
of  Eqns  (F-6)  to  (F-10),  (F-31)  to  (F-34),  (F-43),  and  Eqn  (9)  with  the 
dielectrics  from  Eqn  (20),  (32),  or  (F-52)  are  solved  by  a set  of 
computer  routines  controlled  by  a master  program,  IONION.  IONION  sets 
up  the  problem  by  reading  in  the  initial  distribution  parameters  and 
initializing  the  variables  needed  by  the  solution  routines.  It  then 
controls  the  solution  of  the  problem  by  calling  the  proper  routines  in 
turn  finally  terminating  when  the  solution  is  complete.  The  program  is 
also  set  up  to  handle  magnetic- field- free  electron-ion  instabilities  such 
as  the  ion  acoustic  instability,  but  this  feature  has  never  been  tested. 

The  unmagnetized  versus  the  magnetized  approach  is  selected  by 
setting  the  input  magnetic  field  to  zero  or  not.  Once  this  selection  is 
made,  the  program  automatically  uses  the  proper  form  of  the  electron 
dielectric. 

A block  diagram  of  the  entire  program  showing  the  routines  called 
by  IONION  in  the  order  they  are  called  is  in  Fig.  G-l.  The  figure  shows 
all  of  the  routines  or  groups  of  routines  used  in  the  solution.  Each  of 
these  routines  is  discussed  below  in  the  order  in  which  it  is  first 
called.  The  discussions  Include  brief  explanations  on  what  the  routines 
do  and  how  they  work. 

The  first  routine  called  is  PLASMA.  This  routine  sets  the  reference 
frame  and  calculates  the  fundamental  plasma  parameters  needed  throughout 
the  rest  of  the  program.  The  reference  frame  is  that  shown  in  Fig.  6. 

The  next  routine  is  GRIDLOC  which  locates  the  region  in  the  first 
quadrant  of  the  kx-k2  plane  where  significant  wave  growth  occurs.  It 
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attempts  to  locate  in  k^,  k^,  and  in  the  radial  direction  as  precisely 
as  possible  the  exact  boundaries  of  the  significant  growth  region. 

GRIDLOC  calls  MAXROCT  which  controls  the  linear  dispersion  relation 
solution  for  each  trial  kx  and  k2  from  GRIDLOC.  Of  the  wave  modes  found, 
MAXROOT  saves  the  one  with  the  largest  growth  rate.  Usually,  the  Input 
parameters  are  adjusted  so  that  the  root  finding  routines  look  only  for 
growing  modes. 

The  routine  SETUP  is  called  by  MAXROOT  to  convert  parameters  such 
as  wave  k and  Debye  length  from  rectangular  to  polar  form.  This  form 
Is  needed  by  the  dispersion  relation  routines  called  by  the  root  finders. 

The  next  set  of  routines,  CONTOUR,  SEEK,  and  ANGLER,  called  by 
MAXROOT  locates  the  dispersion  relation  roots  by  searching  over  a large 
area  of  the  complex  frequency  plane.  The  set,  developed  by  Girl  and 
Baum  (Ref  24),  divides  up  the  plane  into  a series  of  boxes,  does  a 
contour  Integration  around  each  box,  and  then  uses  an  extension  of  the 
Cauchy  residue  theorem  to  locate  the  roots  in  each  box. 

MAXROOT  then  calls  CROOT  which  uses  a second  order  Newton- Raphson 
Iteration  scheme  to  improve  on  the  roots  located  by  the  CONTOUR  set. 

CROOT  needs  reasonably  accurate  initial  guesses  so  it  cannot  be  used 
alone.  CROOT  is  also  called  later  In  the  program  to  periodically  update 
the  wave  modes  of  the  evolved  distributions.  In  this  capacity.  It  uses 
the  previously  calculated  modes  as  Initial  guesses. 

The  root  finding  routines  call  CFCTS  which  is  the  analytical  disper- 
sion relation.  CFCTS  uses  the  analytical  form  of  the  ion  and  electron 
dielectrics  given  by  Eqns  (20)  and  (32).  Therefore,  CFCTS  can  only  be 
used  in  the  location  of  the  roots  of  the  initial  Maxwellian  distributions. 

CFCTS  calls  the  set  of  routines  DISP,  DISPPR,  DISPA,  DISPI,  DISPF, 
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DISPS,  GQINTM,  and  ERF  which  evaluates  the  plasma  dispersion  function. 
This  is  the  function  tabulated  by  Fried  and  Conte. 

If  the  magnetized  approach  has  been  selected,  CFCTS  calls  ELDIEL 
which  computes  the  magnetized  version  of  the  electron  dielectric,  Eqn 
(20).  ELDIEL  also  calls  the  DISP  package. 

When  GRIDLOC  Is  finished  having  determined  the  boundaries  of  the 
k plane  grid,  control  returns  to  the  main  program,  IONION.  IONION  then 
calls  the  routine  GRIDS  which  calculates  the  k plane  grid  for  the  diffu- 
sion coefficient  integrals  and  the  velocity  mesh  for  the  species  distri- 
butions. All  three  distributions  use  the  same  velocity  mesh,  and  their 
corresponding  transformations  are  described  In  Appendix  F.  Thus,  the 
same  mesh  point  refers  to  three  different  velocities  depending  on  which 
distribution  Is  being  referenced.  The  k grid  and  velocity  mesh  are 
uniform  and  held  fixed  throughout  the  rest  of  the  program. 

BOSS  is  the  next  routine  called.  This  routine  times  an  activity 
and  is  called  once  when  entering  the  activity  and  again  when  leaving  it. 
The  routine  determines  the  amount  computer  time  remaining  when  entering 
the  activity  and  terminates  the  program  execution  If  the  activity  is 
expected  to  take  too  long.  For  this  determination,  BOSS  keeps  a running 
average  of  the  time  the  activity  takes.  In  IONION,  the  activity  timed 
Is  one  phase  of  the  evolution. 

Next,  MAXROOT  Is  called  directly  from  IONION  for  each  wave  k on  the 
grid.  The  fastest  growing  mode  for  each  k Is  saved  for  the  evolution 
calculation.  If  no  suitable  mode  is  found  at  a given  k,  a dummy  root 
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Is  substituted.  The  dumnjy  roots  zero  the  values  of  the  diffusion  coeffl- 

.Jk 

clent  Integrands  at  their  corresponding  locations  on  the  k grid.  Hence, 
only  significantly  growing  modes  will  build  up  In  energy  and  cause  diffu- 
sion of  particles.  This  part  of  the  program  completes  the  linear  Insta- 
bility problem  solution. 

■ 

The  next  routine  called  Is  BOLTZMN  which  places  the  initial  Maxwell- 
Boltzman  species  distributions  In  their  arrays.  Each  value  of  a distri- 
bution In  Its  array  corresponds  to  a particular  velocity  mesh  point. 

VELLOC  is  called  next  to  locate  the  evolution  control  velocity  mesh 
point.  This  point  Is  two  thermal  widths  from  the  initial  maximum  of  the 
air  distribution  and  is  at  the  same  angle  to  the  negative  vx  axis  as  the 
propagation  angle  of  the  initially  fastest  growing  wave  mode.  The  dis- 
tribution slope  at  this  point,  which  lies  In  the  upper  resonant  region 
shown  in  Fig.  4,  Is  used  to  control  the  evolution.  The  less  dense  air 
distribution  Is  used  becau-3  It  evolves  faster.  FSLOPE  Is  then  called 
to  compute  or  depending  on  whether  the  propagation  angle  Is 

greater  than  45  degrees  or  not.  The  change  In  the  distribution  slope  Is 
a good  measure  of  the  change  In  the  growth  rate  of  the  fastest  growing 
mode.  The  slope  computed  at  this  time  Is  periodically  compared  with  the 

slope  during  the  evolution  of  the  air  distribution.  When  the  slope 

l 

change  becomes  too  large,  the  air  evolution  Is  stopped  so  the  debris 
evolution  can  be  brought  up  to  date  and  the  wave  modes  can  be  updated. 

The  restart  file  routines  SAVEIT,  ALTER,  and  EXTEND  are  called  next. 

1 

The  program  can  be  restarted  after  the  linear  solution  or  after  any 
phase  of  the  evolution  using  the  Information  on  the  restart  file.  The 
* plotting  routine,  FPLOT,  Is  also  called  at  this  time.  This  routine  calls 

several  plotting  routines  In  the  AFWL  computer  system  library,  METALIB. 

I ; 

■; 
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The  resulting  plots  are  three-dimensional  perspective  views  of  the  plasma 
distributions. 

CONSERV  is  called  next.  This  routine  checks  for  conservation  of 
matter,  momentum,  and  energy.  It  also  calculates  the  individual  velocity 
moments  of  the  distributions.  It  is  called  in  each  phase  of  the  evolu- 
tion to  check  on  the  progress  of  the  program. 

The  next  routine,  EVOLVE,  controls  the  solutions  of  the  species  dif- 
fusion equations.  These  solutions  are  the  time  evolutions  of  the  species 
distributions.  EVOLVE  first  calls  ARRAY  for  the  air  distribution  to 


' t 


convert  its  two-dimensional  format  into  the  one-dimensional  format  needed 
by  the  method  of  lines.  It  then  calls  up  the  GEARB  package  by  calling 
DRIVEB  and  controls  the  progress  of  the  air  evolution  by  periodic  calls 
to  FSLOPE.  When  the  current  phase  of  the  air  evolution  Is  finished, 

ARRAY  is  called  again  to  put  it  back  Into  its  original  two-dimensional 
format.  Then,  the  debris  distribution  Is  evolved  In  a similar  manner  up 
to  the  same  ending  time  as  the  air  by  successive  calls  to  ARRAY,  DRIVEB, 
and  ARRAY.  This  step  completes  a phase  of  the  evolution. 

The  GEARB  package  consists  of  the  routines  DRIVEB,  INTERP,  STIFFB, 
COSET,  PSETB,  DECBR,  and  SOLBR.  These  routines  solve  the  ordinary 
differential  equation  set  derived  in  Appendix  F.  These  routines  in  turn 
call  the  set  of  routines  DIFFUN,  PDB,  COEF,  DIFCOEF,  and  COEFINT  which 
define  the  equation,  the  coefficients,  and  the  Jacobian.  The  Inte- 
grands for  the  diffusion  coefficients  are  constructed  in  DIFCOEF,  and  the 
integrations  are  done  by  COEFINT.  COEFINT  Is  designed  to  run  extremely 
rapidly  because  millions  of  diffusion  coefficient  integrals  must  be  done 
in  a typical  evolution. 

When  a phase  of  the  evolution  is  complete;  CONSERV,  SAVEIT,  FPLOT, 
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and  BOSS  are  called  to  check,  save,  and  output  the  results.  BOSS  Is 
called  again  to  start  timing  the  next  phase. 

The  phase  begins  by  the  recalculation  of  the  wave  modes  using  NXTRTS. 
Because  the  ion  distributions  are  no  longer  Maxwellian,  the  numerical 
form  of  their  dielectrics  must  be  used.  NXTRTS  controls  this  updating 
by  calls  to  CROOT  which  locates  the  new  roots.  CROOT  calls  the  routine 
GENDR  which  is  the  dispersion  relation.  GENDR  calls  DRTERM  which  does  the 
Ion  dielectric  integrals,  Eqn  (F-50).  The  electron  dielectric  is  com- 
puted by  GENDR  using  the  same  procedures  which  are  used  by  CFCTS. 

FSLOPE  is  then  called  to  update  the  slope  of  the  air  distribution. 
Finally,  the  anomalous  transport  coefficients  are  computed  using  ANOMAL. 

If  significantly  growing  wave  modes  still  exist  and  the  transport  coef- 
ficients are  still  increasing,  the  program  loops  back  to  the  EVOLVE  call 
to  continue  this  phase  of  the  evolution.  Otherwise  execution  terminates 
as  the  program  is  finished. 

Two  groups  of  utility  routines  are  called  throughout  the  program. 

The  first  group  is  the  SNAP  set  which  outputs  the  value  of  a variable  at 
the  time  the  call  to  SNAP  Is  made.  This  set  is  a very  useful  diagnostic 
aid.  The  other  group  Is  SINGRT,  SINTRP,  and  XYZSPLN.  These  routines 
fit  a smooth  cubic  polynomial  surface  to  a function  defined  by  values  In 
an  array.  The  routines  then  Integrate,  differentiate,  or  Interpolate  on 
that  surface.  All  of  the  program  numerical  Integrations  except  those  for 
the  diffusion  coefficients  are  done  in  this  manner. 
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Appendix  H 

Tables  of  the  Dissertation  Plasma  Distributions 

Below  are  tabulated  all  of  the  plasma  distributions  used  In  the 
dissertation.  The  parameters  listed  In  the  tables  are  as  follows: 
n - number  density 
q - units  of  charge 

AMU  - atomic  mass  units  (atomic  weight) 

V - velocity  (In  the  computer  reference  frame  shown  in  Fig.  6) 
Tx  - temperature  In  the  x direction 
Tz  - temperature  In  the  z direction 


'Dx 


Pz 


Debye  length  In  the  x direction 
Debye  length  in  the  z direction 


top  - plasma  frequency 
_n  - gyro  frequency 


CaJ 


LH  - lower  hybrid  frequency 
- Larmor  radius  (gyro  radius) 
ion  acoustic  speed 


ac 


- Alfven  speed 

cx  - thermal  speed  In  the  x direction 
cz  - thermal  speed  In  the  z direction 
VR  - relative  velocity  between  the  debris  and  air 
B - magnetic  field 
^ - electron  beta 

A given  distribution  is  uniquely  determined  by  the  selection  of  densi- 
ties, charges,  atomic  weights,  and  temperatures  for  the  debris  and  air 
by  the  selection  of  the  electron  temperatures;  and  by  the  selection  of 
the  relative  velocity  and  magnetic  field  strength.  The  remaining  para 
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meters  are  computed  assuming  charge  neutrality,  zero  current,  and  the 
Individual  debris,  air,  and  electron  distributions  are  Maxwellian. 

The  air  densities  and  air  atomic  weight  In  the  last  four  tables. 
Tables  H-VI  to  H-IX,  are  representative  of  an  assumed  predominately 
atomic  oxygen  atmosphere  of  the  earth  existing  at  200  and  600  km.  The 
debris  atomic  weight  In  these  four  tables  Is  that  of  aluminum  which  con- 
veniently represents  an  average  of  all  the  materials  In  the  debris  of  a 
high  altitude  nuclear  explosion.  The  remaining  parameters  which  uniquely 
determine  these  simulation  distributions  are  chosen  to  approximate  the 
conditions  which  could  be  found  In  the  debris-air  mixing  region  shown  in 
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Table  H-I 

Equal  and  Opposite  Cold  Beam  Distribution 


n 

q 

AMU 

V 

(cm-3) 

(units) 

(AMU) 

(cm/sec) 

Debris 

5.1(1010) 

1 

20 

-9.9(106) 

Air 

5.0(1010) 

1 

20 

N 

O 

H 

O 

• 

H 

Electron 

i.o(ion) 

- 

- 

0 

Tx 

Tz 

*DX 

** 

<eV) 

(eV) 

(cm) 

(cm) 

Debrl s 

1 

1 

3. 3(10~3) 

3. 3(10'3) 

Air 

1 

1 

3.3(1 0—  3 ) 

3.3(10'3) 

Electron 

2.8(104) 

2.8(104) 

3.9U0"1) 

3. quo"1) 

% 

_n_ 

^LH 

J°L 

(sec'1) 

(sec"1) 

(sec-1) 

(cm) 

Debris 

6.7 (107) 

0 

0 

OO 

Air 

6 . 6 { 107) 

0 

0 

oO 

Electron 

1.8(1010) 

0 

- 

oo 

vac 

valf 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

2 . 6 (107 ) 

0 

2.2(105) 

2.2 (105 ) 

Air 

2 . 6 ( 107 ) 

0 

2.2(105) 

2.2(105) 

Electron 

- 

- 

7.0(109) 

7 . 0 ( 109 ) 

VR  (cm/sec)  = 2.0(10^) 
B (guass)  = 0 


/e  ’ 1:0 
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Table  H-II 

Bump-on-the-Tail  Distribution 


t 


n q AMU  V 

(era-3)  (units)  (AMU)  (cm/sec) 

Debris  1.0(10n)  1 20  -5.0(105) 

Air  1.0(109)  1 20  5.0(107) 

Electron  1.0(10n)  - n 0 


(eV) 


Debris 

103 

103 

7 . 4 ( 10— 2 ) 

7.4(  0-2) 

Air 

102 

102 

2.4(10-1) 

2.4(10-1) 

El ectron 

2.8(104) 

2.8(104) 

CO 

• 

KD 

1— » 
O 

1 

I— * 

3.9(10-1) 

~P 

_n. 

~LH 

A. 

(sec-1) 

(sec-1) 

(sec-1) 

(cm) 

Debris 

9.3(107) 

0 

0 

-e> 

Air 

9.3(106) 

0 

0 

c*> 

Electron 

1.8(1010) 

0 

- 

vac 

valf 

Cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

3.9(107) 

0 

6.9(106) 

6.9(106) 

Air 

5.3(106) 

0 

2.2 ( 106 ) 

2 . 2 ( 106 ) 

El ectron 

- 

- 

7.0(109) 

7.0(109) 

VR  (cm/sec)  = 5.0(107) 


B (guass)  * 0 

fit  * “ 


4 


I 

i 
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Table  H-III 

Magnetized  Equal  Beam  Distribution 


n 

q 

AMU 

V 

(cm-3) 

(units) 

(AKU) 

(cm/sec) 

Debris 

5.1(1010) 

1 

20 

-9.9(106) 

Air 

5.0(1010) 

1 

20 

1.0(107) 

Electron 

1.0(1011) 

- 

- 

0 

Tx 

Tz 

*Dx 

>Dz 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

1 

1 

3.3(10"3) 

3.3(10"^) 

Air 

1 

1 

3.3(10'3) 

3.3(10"3) 

Electron 

10"  2 

10"2 

2.3(10'4) 

2.3(10'4) 

% 

-TL 

^LH 

J°L 

(sec'1) 

(sec"1) 

(sec-1) 

(cm) 

Debris 

6.7(107) 

4.8(105) 

4.7(107) 

2.K101) 

Air 

6.6 { 107) 

4.8(105) 

4.7(107) 

2.K101) 

Electron 

1.8(1010) 

1.8(1010) 

- 

3.4(10”4) 

vac 

valf 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

3.8(105) 

2 . 2 ( 108 ) 

2.2(105) 

2.2(105) 

Air 

3.8(105) 

2 . 2 ( 108 ) 

2.2(105) 

2.2(105) 

Electron 

- 

- 

4.2(106) 

4.2(106) 

VR  (cm/sec)  = 2.0(107) 
B (guass)  = 1.0(103) 


* 4.1(10-8) 
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Table  H-IV 


Initial 

Bump-on-the-Tail 

Distribution 

for  the  Test 

Evolution 

n 

q 

AMU 

V 

(cm-3) 

(units) 

(AMU) 

(cm/sec) 

Debri s 

1.0(H)1*, 

1 

20 

-3.0(10S) 

Air 

1.0(109) 

1 

20 

3. 0 ( 107) 

Electron 

1.0(H)11) 

- 

- 

0 

Tx 

Tz 

*Dx 

*Dz 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

103 

103 

7.4(10"2) 

7.4(10"2) 

Air 

102 

102 

2 . 4 ( 10“ 1 ) 

2.4(10"1) 

Electron 

2.8(104) 

2.8(104) 

3.9(H)-1) 

3.9(H)-1) 

% 

_n_ 

^LH 

Al 

(sec-1) 

(sec"1) 

(sec"1) 

(cm) 

Debris 

9.3(107) 

0 

0 

06 

Air 

9.3(106) 

0 

0 

oO 

Electron 

1 . 8 { 1010 ) 

0 

- 

06 

vac 

valf 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

3.9(107) 

n 

6.9(106) 

6 . 9 ( 106 ) 

Air 

5.3(106) 

0 

2.2(10®) 

2 . 2 ( 106 ) 

Electron 

- 

- 

7 . 0 ( 109 ) 

7 .0 ( 109) 

VR  (cm/sec)  = 3.0(10^) 
B (guass)  = 0 


A 
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Table  H-Y 

Large  Relative  Velocity  Bump-on-the-Tatl  Distribution 


n 

q 

AMU 

V 

(cm"3) 

(units) 

(AMU) 

(cm/sec) 

Debris 

i.o(ion) 

1 

20 

-9.9(105) 

Air 

1.0(109) 

1 

20 

9 . 9 ( 107) 

Electron 

l.Kio11) 

- 

- 

0 

Tx 

T2 

*0* 

*DZ 

(eV) 

(eY) 

(cm) 

(cm) 

Debris 

10* 

103 

7.4(10"2) 

CM 

1 

O 

r—i 

** 

• 

is 

Air 

CNJ 

o 

102 

2.4(10'1) 

2. 4( 10"1 ) 

Electron 

2.8(104) 

2.8(104) 

3.  gfio'1) 

s.gdo"1) 

% 

60  LH 

J\ 

(sec"1) 

(sec'1) 

(sec-1) 

(cm) 

Debris 

9.3(107) 

0 

0 

Air 

9. 3( 106) 

0 

0 

C*> 

Electron 

1.8(1010) 

0 

- 

OO 

vac 

valf 

cx 

C2 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

3.9(107) 

0 

6.9(10®) 

6.9(10®) 

Air 

5.3(10®) 

0 

2 .2 ( 106) 

2. 2 ( 106) 

El ectron 

- 

- 

7.0(109) 

7.0(109) 

(cm/sec)  a 1.0(10®) 


B (guass)  * 0 

A ■ “ 
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Table  H-VI 

Low  Relative  Velocity  Distribution  Simulating  200  km 


n 

q 

AMU 

V 

(cm"3) 

(units) 

(AMU) 

(cm/sec) 

Debris 

2.0(1016) 

1 

27 

-6.7(106) 

Air 

1.0(1010) 

1 

16 

1.3(107) 

Electron 

3.0(1010) 

- 

- 

0 

T 

X 

Tz 

^Dx 

*Dz 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

5.0(10*) 

5.0(102) 

1 . 2 ( 10~ 1 ) 

1.2(H)'1) 

Air 

5.0(102) 

5 . 0 ( 102 ) 

1.7(10_1) 

1 . 7 ( 10~ 1 ) 

Electron 

1.0(104) 

1.0(104) 

4 . 3( 10— 1 ) 

4 . **  * 10" 1 ) 

A 

_TL 

Ah 

A 

(sec"1) 

(sec-1) 

(sec-1) 

(cm) 

Debris 

3.6(107) 

0 

0 

oO 

Air 

3.3(107) 

0 

0 

o* 

Electron 

9.8(109) 

0 

- 

CO 

vac 

val  f 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

1.7(107) 

0 

4.2(106) 

4.2(106) 

Air 

1 . 7 ( 107 ) 

0 

5 . 5 ( 106 ) 

5.5(106) 

Electron 

- 

- 

4.2(109) 

4.2(109) 

VR  (cm/sec)  = 2.0(107) 
B (guass)  = 


0 
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Table  H-V1I 

High  Relative  Velocity  Distribution  Simulating  200  km 


n 

q 

AMU 

V 

(cm-3) 

(units) 

(AMU) 

(cm/sec) 

Debrl s 

2.0(1010) 

1 

27 

-2.3(107) 

Air 

1.0(1010) 

1 

16 

4. 7(107) 

El ectron 

3.0(1010) 

- 

- 

0 

T 

X 

Tz 

*Dx 

*Dz 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

1.0(103) 

1.0(103) 

1.7(10"*) 

1 . 7 ( 10“ 1 ) 

Air 

1.0(103) 

1.0(103) 

2 . 4 ( 10” 1 ) 

2.4(10"*) 

El ectron 

1.0(104) 

1.0(104) 

4.3(10"*) 

4.3(10"*) 

% 

—ft. 

(sec-*) 

(sec-1) 

(sec-1) 

(cm) 

Debris 

3.6(107) 

0 

0 

oO 

Air 

3.3(107) 

0 

0 

DO 

Electron 

9.8(109) 

0 

- 

oo 

vac 

Valf 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

1. 9(107) 

0 

6.0(106) 

6. 0 ( 106 ) 

Air 

2 . 0 ( 107 ) 

0 

7 . 8 ( 106 ) 

7 . 8 ( 106 ) 

Electron 

- 

- 

4 . 2 ( 109 ) 

4.2(109) 

VR  (cm/sec)  = 7.0(107) 


B (guass)  = 0 

fit  - ~ 

L 
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Table  H-VIII 

Low  Relative  Velocity  Distribution  Simulating  600  km 


n 

q 

AMU 

V 

(cm-3) 

(units) 

(AMU) 

(cm/sec) 

Debris 

2 . 0( 107) 

1 

27 

-6 . 7 ( 106 ) 

Air 

1.0(107) 

1 

16 

1 . 3 ( 107 ) 

Electron 

3 . 0 ( 107 ) 

- 

- 

0 

Tx 

Tz 

>D* 

^DZ 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

5.0(102) 

5 . 0 ( 102 ) 

3.7 

3.7 

Air 

5.0(102) 

5 . 0 ( 102 ) 

5.3 

5.3 

El ectron 

1.0(1C3) 

1.0(104) 

4.3 

1.4(H)1) 

A) 

P 

_n 

^LH 

(sec-1) 

(sec-1) 

(sec-1) 

(cm) 

Debris 

1.1(106) 

1.8(104) 

1.1(106) 

3. 7 ( 102 ) 

Air 

1.0(106) 

3.0(104) 

9.9(105) 

4.4(102) 

Electron 

3.1(108) 

8.8(108) 

- 

2.1 

vac 

val  f 

cx 

cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cm/sec) 

Debris 

1 . 4 ( 107 ) 

4.7(108) 

4 . 2 ( 106 ) 

4.2(106) 

Air 

1.4(107) 

8 . 7 ( 108) 

5 . 5 ( 106 ) 

5.5(106) 

Electron 

- 

- 

1.3(109) 

4.2(109) 

VR  (cm/sec)  = 2.0(107) 
B (guass)  = 5.0(H)1) 
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Table  H-1X 

High  Relative  Velocity  Distribution  Simulating  600  km 


n 

q 

AMU 

V 

(cm“3) 

(units) 

(AMU) 

(cm/sec) 

Debris 

2.0(107) 

1 

27 

—2 . 3 ( 107 ) 

Air 

1.0(107) 

1 

16 

4 . 7 ( 107 ) 

Electron 

3. 0 ( 107) 

- 

- 

0 

Tx 

Tz 

^Dx 

^Dz 

(eV) 

(eV) 

(cm) 

(cm) 

Debris 

1.0(103) 

1.0(103) 

5.3 

5.3 

Air 

1.0(103) 

1.0(103) 

7.4 

7.4 

Electron 

1.0(103) 

i.o(io4) 

4.3 

l.AdO1) 

% 

^°LH 

Si 

(sec-1) 

(sec"1) 

(sec-1) 

(cm) 

Debris 

l.Kio6) 

1.8(104) 

1.1(106) 

1.3(103) 

Air 

1.0(106) 

3. 0 ( 104 ) 

9.9(105) 

1.5(103) 

El ectron 

3. 1(108) 

8.8(108) 

- 

2.1 

Vac 

valf 

cx 

Cz 

(cm/sec) 

(cm/sec) 

(cm/sec) 

(cir./sec) 

Debris 

1.5(107) 

4.7(108) 

6.0(106) 

6.0(106) 

Air 

1.7 (107) 

8.7(108) 

7 . 8 ( 106 ) 

7.8(106) 

Electron 

- 

- 

1.3(109) 

4.2(109) 

VR  (cm/sec)  = 7.0(107) 
B (guass)  c 5.0(10^) 
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Appendix  I 

Evolution  of  the  Bump-on-the-Tail  Instability 


The  analytical  quasilinear  solution  for  the  bump-on- the-tall 
Instability  Is  done  using  a combined  debris  and  air  Ion  distribution. 

The  same  procedures  are  used  which  are  employed  by  Davidson  (Ref  14:  174- 
182)  In  his  solution  of  the  equivalent  electron  bump-on-the-tall  problem. 
For  the  test  Ion  problem  shown  In  Table  H-IV,  the  electron  temperature  is 

O A 

large  enough  so  the  electron  dielectric  which  goes  like  1/k (See 
Eqns  (53),  (57),  and  (58).)  Is  less  than  one.  Therefore,  the  electrons 
can  be  ignored  because  their  contribution  to  the  dispersion  relation, 
and  hence  to  the  unstable  waves  causing  the  time  evolution  of  the  ion 
distribution,  is  negligible. 

The  equations  of  interest  for  the  ions  are  from  Eqns  (10),  (11), 
and  (12): 


- U 3-^]  (M) 


3 T i ^ (A  fGA  ^ 


(1-2) 


r\/  ,/N  O -r>  3^  f //  £ ^ ^ ^ - (1  31 

= nr 

where  the  subscript  o<  has  been  dropped  and  the  equations  are  written  in 
a one-dimensional  form. 

The  diffusion  coefficient  has  two  different  approximate  forms  with- 
in and  outside  the  resonant  region.  In  the  resonant  region, 


Jisv C'bj  ^ ^ 


(1-4) 


so  the  corresponding  coefficient  is 
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x. 

)M  Z(AJ)  S(J>*r-  c^(A,y))  (I- 


5) 


where  Eqn  (F-37),  60^ (k,t)  t 0,  and 

fg£  


V-SCJur-ctfA^)  (1-6) 


have  been  used.  Outside  the  resonant  region, 

[ Jlyxr  — a<*.  (Ut-jjf}  ] > ">  £*>U  (>4, 

so  the  nonresonant  form  of  the  coefficient  is 

d -'i  - is-^Tjjl  UAJli^LLil  _ 


(1-7) 


(1-8) 


The  resonant  form  after  the  symmetries  in  Fqns  (A-24)  and  (A-26)  are 
applied  is 


Dsx  \ cUt.  jO<rC-^v — ^ ( -4, (1-9) 

o 

This  integration  is  done  by  assuming 

^ c^po  (I- 10) 

s°  [ r ca,/)\  j j?  o d-ii) 

“4  ■ 

This  diffusion  coefficient  is  negligible  outside  the  resonant  region 
because  the  wave  spectral  energy  does  not  grow  above  the  background 
thermal  levels. 

The  approximation  in  Eqn  (I- 10)  is  valid  for  long  wavelength  modes 
so  all  but  the  first  term  on  the  right  hand  side  of 

~ - ( I + ^ j*-  + " ' ^ (1-12) 

can  be  Ignored.  This  is  the  real  part  of  the  Bohm-Gross  dispersion 
relation  which  is  derived  for  electrons  in  Krall  and  Trivelpiece  (Ref 
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19:  383-392).  The  corresponding  Imaginary  part  Is 


Tr  . . c)/e(sVijf') 


(1-13) 


The  long  wavelength  approximation  Is  a good  one  for  the  test  distribution 


as 


Da 


AX 


S'  ( lo~x>)  <<  I 


(1-14) 


for  the  growing  modes  of  Interest. 

The  resonant  form  of  Eqn  (1-1)  Is 


J I (1-15) 

where  Dr  Is  given  by  Eqn  (I- 11) - The  time  derivative  of  Eqn  (1-11) 
combined  with  Eqns  (1-2),  (1-10),  and  (1-13)  Is 

, -Y  f . t T 1«-\ 

■jj — = 1 r^ro^r  t-^F  J (1-16) 


aft, 


where  Eqn  (1-11)  has  been  used  again  to  eliminate  the  spectral  energy 
term.  The  integral  of  Eqn  (1-16)  Is  j. 

~ - o')  JL 

When  Eqns  (1-15)  and  (1-16)  are  combined  and  integrated. 


(1-17) 


= a [ (Ms) 


c)y\r  L Cu’fo'V-  d-f 


(1-19) 


/ / _ _ / s y.  \ , ^ r O/t.  Cyv',  ^ ~ «y~} 

J 

Although  Eqns  (1-17)  and  (1-19)  are  the  resonant  solution,  they  still 
are  not  in  a form  suitable  for  a hand  calculation.  The  best  approach  is 
to  find  the  time  asymptotic  solution. 
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The  first  step  In  the  asymptotic  solution  Is  to  multiply  Eqn  (1-15) 
by  f0(v,t)  and  Integrate  with  respect  to  v: 

& |<£r  [ /.  c-,»f = - d^c^I  ^~\ - o (i-2i) 


a 

-<*»  ~o» 

The  Integrated  part  resulting  from  the  right  hand  side  of  Eqn  (1-20) 
vanishes  because  resonant  diffusion  is  negligible  outside  the  resonant 
region.  Because  the  time  rate  of  change  of  a positive  quantity  Is 
negative  for  all  time, 

„ . , . r ( yv-  jt  _ > oo'y  1 t 

£r J = ° (1-22) 


Either  or  both  terms  can  be  zero.  If 


then  from  Eqn  (1-17), 


O^C^j  J - 


(1-23) 


(1-24) 


so  the  growth  rate  from  Eqn  (1-13)  becomes  negative.  Equations  (1-19) 
and  (1-23)  also  Imply 

/•o,  -*-»  - /« (I'25) 

Because  the  particle  energy  is  much  larger  than  the  Initial  wave  energy, 
the  Initial  diffusion  should  be  small.  Therefore,  the  second  term  can 
be  neglected  and 

<V.  0 (1-26) 


1 e:e. 


(1-26) 
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Equations  (1-24)  and  (1-26)  contradict  so  Eqn  (1-23)  must  be  Incorrect 


— o 


(1-27) 


Thus,  the  asymptotic  behavior  of  the  distribution  In  the  resonant  region 
Is  to  form  a plateau  with  the  wave  growth  rates  going  to  zero. 

The  nonresonant  diffusion  coefficient,  Eqn  (1-8),  is  used  to  deter- 
mine the  behavior  of  the  bulk  of  the  distribution  (the  debris  distribu- 
tion). Because  the  phase  velocities  of  the  growing  waves  lie  In  the 
resonant  velocity  range  which  Is  much  greater  than  the  thermal  velocity 


of  the  bulk  distribution. 


[to 


(1-28) 


Therefore,  Equation  (1-8)  becomes 


d-29) 


bpH.fajS)  - Jjt  \ <Lk 


(1-30) 


Note  the  coefficient  Is  Independent  of  velocity.  The  nonresonant  form 


of  Eqn  (1-1)  Is  then: 


(1-31) 


When  T(t)  Is  defined  as 


T(jf)  - 


(1-32) 


Equation  (1-31)  becomes 


(1-33) 
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The  Green's  function  solution  of  this  equation  with  an  Initial  Maxwellian 
distribution  for  the  bulk  plasma  Is 


!/  C s>v\sv~t'  ^ 

y / S. HI } *■  ~ (/2-LxT  (1-34) 

■fj'*',*)  ^ l vn-LxT  + lv>:>3)  JL  K } 

where  T Is  the  bulk  Ion  temperature.  Note  that 

roo  >>  TCo)  -Tor  st  •>>  c*Jj  (1-35) 

has  been  assumed.  The  nonresonant  solution,  Eqns  (1-32)  and  (1-34), 
shows  that  the  effective  temperature  of  the  bulk  of  the  Ions  Increases 
as  the  wave  energy  Increases.  This  Increase  Is  due  to  bulk  plasma 
oscillations  Induced  by  the  electrostatic  waves.  A sketch  of  the  final 
time  asymptotic  Ion  distribution  showing  both  the  resonant  and  non- 
resonant changes  Is  Included  In  Fig.  7. 

The  growth  of  the  wave  energies  and  the  growth  of  the  bulk  plasma 

oscillations  Is  fed  by  a loss  in  the  energy  of  the  resonant  particles 

in  the  bump.  In  fact,  (Ref  19:  527-532) 

change  in  change  change  In 

nonresonant  _ In  total  L y resonant  /T 

particle  ~ wave  ■%,  A particle  ' ’ 

energy  energy  energy 

‘ i • 


I 
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Appendix  J 

Tables  of  the  Linear  Dispersion  Relation  Solutions 

The  linear  dispersion  relation  solutions  for  growing  wave  modes 
are  tabulated  for  all  of  the  distributions  which  are  evolved.  The 
solutions  are  for  Maxwellian  debris,  air,  and  electron  distributions  so 
they  are  the  wave  modes  which  exist  at  the  beginning  of  an  evolution. 

An  additional  dispersion  relation  solution  table.  Table  J-V,  is  included 
to  show  a possible  mixture  of  modified  two-stream  and  ion-ion  modes. 

This  table  is  discussed  In  the  section  on  magnetized  evolutions  in 
Chapter  IV. 

In  all  of  the  tables,  only  the  growing  wave  modes  in  the  first 


quadrant  of  the  kx-k£  plane  where 

0 < K 

(J-l) 

0 — -Jt  j 4 oA 

(J-2) 

are  Included.  The  modes  propagating  In  the  other  quadrants  of  this 
plane  are  related  to  those  tabulated  by  the  equations  In  Appendix  F, 

Eqns  (F-39),  (F-40),  and  (F-42). 

In  the  tables  below,  the  headings  are  the  components  of  the  wave  k 
in  the  x and  z directions;  the  real  and  Imaginary  parts  of  the  wave 
frequency;  and  the  debris,  air,  and  electron  dielectrics  evaluated  at 
the  wave  frequency.  The  dielectrics  are  complex  so  both  the  real  and 
imaginary  parts  are  listed.  The  row  marked  with  an  asterisk  is  the  wave 
mode  with  the  largest  growth  rate.  This  mode  dominates  the  evolution. 


Linear  Solutions  for  the  Bump-on- the-Tall  Test  Case 
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